Comparative Transcriptomics Reveals the Molecular Genetic Basis of Cave Adaptability in Sinocyclocheilus Fish Species

Cavefish evolved a series of distinct survival mechanisms for adaptation to cave habitat. Such mechanisms include loss of eyesight and pigmentation, sensitive sensory organs, unique dietary preferences, and predation behavior. Thus, it is of great interest to understand the mechanisms underlying these adaptability traits of troglobites. The teleost genus Sinocyclocheilus (Cypriniformes: Cyprinidae) is endemic to China and has more than 70 species reported (including over 30 cavefish species). High species diversity and diverse phenotypes make the Sinocyclocheilus as an outstanding model for studying speciation and adaptive evolution. In this study, we conducted a comparative transcriptomics study on the brain tissues of two Sinocyclocheilus species (surface-dwelling species – Sinocyclocheilus malacopterus and semi-cave-dwelling species – Sinocyclocheilus rhinocerous living in the same water body. A total of 425,188,768 clean reads were generated, which contributed to 102,839 Unigenes. Bioinformatic analysis revealed a total of 3,289 differentially expressed genes (DEGs) between two species Comparing to S. malacopterus, 2,598 and 691 DEGs were found to be respectively, down-regulated and up-regulated in S. rhinocerous. Furthermore, it is also found tens of DEGs related to cave adaptability such as insulin secretion regulation (MafA, MafB, MafK, BRSK, and CDK16) and troglomorphic traits formation (CEP290, nmnat1, coasy, and pqbp1) in the cave-dwelling S. rhinocerous. Interestingly, most of the DEGs were found to be down-regulated in cavefish species and this trend of DEGs expression was confirmed through qPCR experiments. This study would provide an appropriate genetic basis for future studies on the formation of troglomorphic traits and adaptability characters of troglobites, and improve our understanding of mechanisms of cave adaptation.


INTRODUCTION
Caves are classified among the special habitat category, and are characterized by extremely scarce of food resources, high concentrations of carbon dioxide, and perpetual darkness. Troglobites exhibit a set of phenotypes and behavioral traits that have arisen from long-term adaptation to these harsh habitat factors (Aspiras et al., 2015). Studies of these traits will contribute not only to improved understanding of the mechanisms involved in their adaptation, but also help to reveal processes underlying speciation, acquisition of novel functions, and the impact of evolution on genetic material.
"Cavefish" mainly refers to fish species that live in subterranean cave water bodies, underground rivers, underground lakes or other similar ecosystems (Borowsky, 2018). Cavefish have developed a series of unique physical and behavioral features in their adaptation to cave habitats, including reduction of eyesight and pigmentation, sensitive sensory organs, unique dietary preferences, and predation behavior (Jeffery, 2009). In recent years, many research related to morphology, molecular biology and ecology has been conducted on various cavefish species such as Astyanax mexicanus (Varatharasan et al., 2009;Hyacinthe et al., 2019) and Phreatichthys andruzzii (Ceinos et al., 2018). However, due to difficulties in adequately sampling wild cavefish species, clear understanding of regulatory mechanisms underlying their adaptive traits requires further investigation (Zagmajster et al., 2010).
As the control center in vertebrates, the brain plays a significant role in regulating body temperature, biorhythm, muscle activity, and the sensory systems. A subset of metazoan genes underlie the biophysical properties and information processing of the brain (Guo, 2003;Drew et al., 2012). The goal of studies of the functions and regulatory mechanisms of these brain genes to understand cell-and circuit-specific effects is ultimately to establish the connections between brain function and adaptive behavior (Croset et al., 2018). For example, Vincent et al. found that the territorial behavior of the round goby (Neogobius melanostomus) is linked to methylation of DNA in the hypothalamus, and the process of DNA methylation synchronized with the appearance of territorial behavior (Somerville et al., 2019). Similarly, by comparing the hypothalamus transcriptomes of domesticated and aggressive undomesticated silver foxes (Vulpes vulpes), Rosenfeld et al. (2019)  which play roles in signaling, development, differentiation, and immunity. These DEGs may be partly responsible for the loss of aggressive behavior after domestication of silver foxes (Rosenfeld et al., 2019).
The teleost genus Sinocyclocheilus (Cypriniformes: Cyprinidae) is a treasured group of wild fish species, endemic to the karst terrain of the east and northwest Yunnan and Guizhou plateau, Guangxi, southwestern China, and encompassing more than 70 known species as of 2019. This genus is the largest known group of cavefishes worldwide, with new species continuing to be reported (Lunghi et al., 2019). Interestingly, the genus Sinocyclocheilus represents an example of two morphotypes (cavefish/surface fish) that have obvious phenotypic and habitat differences existing within one genus. High species diversity with diverse phenotypes makes Sinocyclocheilus as an outstanding model for studying speciation and adaptive evolution, that has attracted much attention in various fields. For example, Yang et al., 2016 conducted comparative genomic and transcriptomic studies on three representative Sinocyclocheilus fish species (surface-dwelling Sinocyclocheilus grahami; semi-cave-dwelling Sinocyclocheilus rhinocerous; and cave-restricted Sinocyclocheilus anshuiensis) and found many genetic changes associated with regressive features, such as gene loss, pseudogenes, mutations, and differential gene expression. The genetic changes identified in this study provide novel insights into cave dwelling adaptations of Sinocyclocheilus fish (Yang et al., 2016). Research into the adaptability of Sinocyclocheilus fishes has largely focused on the most distinct troglomorphic traits (eye degeneration or skin albinism) (Meng et al., 2013), while other organs/tissues such as the brain, muscle, and liver remain insufficiently studied.
In this study, we conducted comparative transcriptomic analysis to examine the transcript profiles from the brains of two sympatric Sinocyclocheilus fish species: surface-dwelling Sinocyclocheilus malacopterus (with complete eye structure) and semi-cave-dwelling S. rhinocerous (with reduced eyesight and pigmentation). The main purposes of this study were to: (i) obtain the high-quality brain transcriptome profiles of cave and surface dwelling Sinocyclocheilus species using Illumina sequencing technology, (ii) identify DEGs by comparing brain transcriptome profiles of cave and surface dwelling species, and (iii) to identify and confirm key candidate genes which related to cave adaptability and development that differ between these two different ecotype fish species through qPCR experiments, and pathway and function annotation of DEGs.
The results of this research are expected to provide valuable information to guide further studies of biocompatibility and improve our understanding of the molecular regulatory mechanisms underlying troglomorphic traits in cavefishes.

Sample Collection
The two fish species S. malacopterus and S. rhinocerous used in this study were captured in Yunnan Province, China, with five and three biological replicates, respectively (Figure 1). The Characteristics of the sampling location are shown in the Table 1.
To avoid RNA degradation, because cavefish die soon after being captured from the cave water body, sterile laboratory appliances were prepared in advance, and all tissue samples were harvested immediately after the fish were collected from the water. Because of the small size of the fish, after general anesthesia with 30 mg/L of MS-222 anesthetic (3-Aminobenzoic acid ethyl ester methanesulfonate) (Sigma-Aldrich, St. Louis, MO, United States), all brain tissues were surgically excised, consolidated into sterile tubes containing RNAlater, and placed in liquid nitrogen. After sacrificing the fish, other tissues (muscles) and organs (heart, liver, skin, etc.) were also collected for other studies. Every effort was made to minimize animal suffering. Samples were stored at −80 • C upon return to the laboratory. All experiments were conducted with approval from the local Ethical committee at Yunnan University in accordance with China's local and global ethical polices (Grant No: ynucae 20190056),  and all procedures were approved and permits granted by the local government.

RNA Extraction, Library Preparation, and Illumina Sequencing
Total RNA was extracted from the brain tissues of each Sinocyclocheilus fish using TRIzol R Reagent (Invitrogen, Carlsbad, CA, United States). A 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, United States) was used to determine RNA purity and integrity, and the RNA samples were quantified using the ND-2000 (NanoDrop Thermo Scientific, Wilmington, DE, United States). High-quality RNA was used to construct the sequencing library. RNA purification, reverse transcription, and library construction were performed according to the manufacturer's instructions (Illumina, San Diego, CA, United States). Brain RNA-seq transcriptome libraries were prepared using an Illumina TruSeq TM RNA sample preparation Kit (San Diego, CA, United States). Poly(A) mRNA was purified from total RNA using oligo-dT-attached magnetic beads, then fragmented in fragmentation buffer. Using these short fragments as templates, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad, CA, United States) with random hexamer primers (Illumina). Synthesized cDNAs were subjected to end-repair, phosphorylation, and "A" base addition according to Illumina's library construction protocol. The libraries were limited to 200-300 bp fragments by gel purification from 2% Low Range Ultra Agarose, followed by PCR amplification using Phusion DNA polymerase (New England Biolabs, Boston, MA, United States) for 15 cycles. After quantification using a TBS380, two RNAseq libraries were sequenced in a single lane on an Illumina HiSeq X ten sequencer (Illumina, San Diego, CA, United States) for 2 × 150 bp paired-end reads. In this study, DNA sequenced using an Illumina HiSeq platform, with assistance from Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China).

De novo Assembly and Sequence Annotation
Clean data from the brain samples were used for de novo assembly using Trinity 1 , after quality control using SeqPrep 2 and Sickle 3 software packages to remove low-quality reads, reads of less than 20 bp in length, and reads containing more than 10% ambiguous bases "N." BLASTX was used to identify those proteins in the NCBI protein non-redundant (NR), and KEGG databases with the highest sequence similarity to each of the assembled transcripts, to gain insight into their function. BLAST2GO 4 (Conesa et al., 2005) was used to obtain GO annotations describing biological processes, molecular functions, and subcellular locations of unique assembled transcripts. Metabolic pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) 5 (Minoru and Susumu, 2000), and functions were predicted using Unigenes by aligning Unigenes sequences to the COG database.

Differential Expression Analysis
To identify DEGs between the two different fish species, the expression level of each transcript was calculated in units of fragments per kilobase of exon per million mapped reads (FRKM) method. RSEM 6 was used to quantify gene and isoform abundances (Dewey and Li, 2011). The R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R), was utilized for differential expression analysis. Genes with expression log2| Foldchange| ≥ 1 and false discovery rates (FDR) ≤ 0.05 were considered to be differentially expressed.

Functional Annotation of Differentially Expressed Genes
Functional-enrichment analyses, including GO and KEGG, were performed to identify which GO terms and metabolic pathways were significantly enriched among identified DEGs, using a Bonferroni-corrected P-value ≤ 0.05 (significantly enriched functional clusters). GO functional enrichment and KEGG pathway analysis were carried out using Goatools 7 and KOBAS 8 , respectively.

Validation of DEGs by Real-Time PCR
In this study, due to quantitative limitations on wild fish brain mRNA samples, we could not validate all of the DEGs that we identified. Therefore, four DEGs were randomly selected from among all target DEGs for validation using real-time PCR analysis. Real-time PCR analysis used the same RNA samples used for transcriptome profiling. Two housekeeping genes (β-actin and GAPDH) were used as internal references. Primer sequences are listed in Table 2, and amplification conditions were as follows: initial denaturation at 95 • C for 5 min, followed by 40 cycles of 95 • C for 5 s, 50 • C (nmnat1)/55 • C (dab1, ophn1, and ubr1) for 30 s, and extension at 72 • C for 40 s. Finally, 60-95 • C was used to produce melting curves for analysis. All reactions were performed with three technical replicates using one biological sample. Relative gene expression was calculated using the 2 − CT method using qPCRsoft3.2 software.

Weighted Gene Co-expression Network Analysis (WGCNA)
Although the quantity of samples in this study did not meet the official recommended standards, in consideration of the rarity of cavefish species and small population sizes, we still launched WGCNA and tried to find additional important genes. We used the R package WGCNA for weighted correlation network analyses. Using gene expression data that assessed transcriptional levels in brain tissue from cave and surface fish, we constructed a co-expression network.
First, we constructed a gene-gene similarity network (Pearson's correlation) for all Unigenes. In this analysis, we filtered Unigenes for those with relative expression levels <1 to calculate optimal power values. We then transformed these values to an adjacency matrix under soft power (beta = 6). Next, all Unigenes were hierarchically clustered and divided into modules. Gene modules corresponding to the branch cutoff of the gene tree were color-coded (networkType: signed; soft power: 6; minModuleSize: 30; minKMEtoStay: 0.3; mergeCutHeight: 0.25). The "gray" module contains Unigenes that cannot be associated with any expression patterns. In this study because we cannot measure physiological data from each wild cavefish species, we have instead produced a simple phenotypic traits table. To find the core genes (hubgenes) that are significantly associated with a phenotypic trait, we focused on the most highly correlated modules (R 2 > | 0.4|, P < 0.005). The most correlated genes, in "blue" and "turquoise, " were selected for further analysis. Information on the functions and related pathways of hub genes were obtained from annotation in GO and KEGG databases.

Sequencing and Assembly
In this study, a total of 425,188,768 clean reads, covering 63,744,868,368 bases were produced form the brain transcriptomes of the two Sinocyclocheilus species (eight specimens) using an Illumina HiSeq 2000 sequencing platform.
Other parameter details of the sequencing data, such as bases number, and GC content, are presented in Table 3. A total of 102,839 Unigenes (unique transcript fragments), with an average length of 1151.19 bp and an N50 of 1734 bp was generated from 184,699 transcripts de novo assembled in Trinity. Information about sequence the length distribution and other details of transcripts and Unigenes is shown in Table 4.

Functional Annotation of Unigenes
Prediction of function, and Unigene categorization were achieved by searching in many databases (  (Figure 2). According to GO terms, 122,719 Unigenes were assigned to GO level 2 terms using GO database annotation under three major functional categories: "biological process, " "cellular component, " and "molecular function." The level 3 GO terms under each functional are shown in Figure 2B. Regulation of biological process (7,646 Unigenes) and organic substance metabolic process (6,517 Unigenes) were the main subcategories under "Biological

DEGs Analysis and Functional Annotation
Sinocyclocheilus malacopterus and S. rhinocerous were designated control group and experimental groups, respectively. A total of 3,289 DEGs were identified, and 2,598 and 691 DEGs, respectively, were found to be respectively, down-regulated and up-regulated in S. rhinocerous (cave-dwelling species). The proportion of down-regulated genes in this surface vs. cave fish comparison. The volcano diagram in Figure 4 provides more detail.
To further filter the candidate DEGs, we attempted to ascertain each DEG's function using GO database annotation. DEGs in the surface-cave fish comparison group were classified into three functional categories in the GO database: biological process, Cellular component, and molecular function. However, several GO terms were over-or under-represented among these DEGs (Figure 5A). In the biological process category, regulation of biological processes and responses to stimuli were overrepresented among the up-regulated gene group, which may likely reflect the highly responsive behavior of cavefish to external stimuli, such as predation. The cells and organelles subcategory of the cellular component category was overrepresented among down-regulated genes. GO functional enrichment analysis of DEGs by using Goatools software, showed that the significant enrichment of DEGs functions was mainly concentrated in two categories: cellular component and molecular function (under which organelle membrane was the most significant subcategory) (P < 0.001) (Figure 5B).

Target DEGs Selection and Expression Trend Validation
After the DEG functional analysis and pathway annotation, 13 of the 3289 total DEGs, putatively related to cave adaptability of the Sinocyclocheilus fish species, were selected. These DEGs are involved in cavefish specific phenotypes and behavioral traits, such as eye, brain, and retina development ( (Table 6). We verified this expression trend through the real-time PCR to ensure the reliability of transcriptome data in this study (Figure 6).

The WCGNA of Cave and Surface Sinocyclocheilus Fish Species
After filtering, 41,958 entered genes were incorporated into 51 color-modules. The module of "turquoise" (12,239 genes) and "blue" (5,989 genes) modules are the two that contain the most genes. Details of the correlation between modules and clustering patterns are shown in Supplementary Material 5. A co-expression network was used to quantitatively link gene expression with two troglomorphic phenotypes, which readily distinguish two types of fish: eyes and horn (Figures 7A,B). In this study, our analysis identified groups of genes correlated with the observed troglomorphic traits. We found that the "turquoise" and "blue" modules showed a significant correlation with differentiation of Sinocyclocheilus fish species according to ecotype (Figures 7C,D). By focusing on these two highly correlated modules, we found tens of hub genes (Figure 8). Functional and pathway annotation indicated that these hub genes with differential expression between surface fish and cavefish are involved in a variety of biological processes, and are divided into the following four major categories: (i) regulation of insulin; (ii) energy and metabolism, especially related to ATP and GTP biosynthesis; (iii) phenotype and behavior; and (iv) nervous system development ( Table 7).

DISCUSSION
Cave-adaptation has ultimately created a unique survival mechanisms and adaptive behaviors in organisms from the chronic stress of harsh habitat factors (Weixian et al., 1997). Many studies have shown that, the complex behaviors that are goal-directed (driven by the "intention" to survive or reproduce, e.g., food hunting, escape from predators) or emotion-related (e.g., aggression, anxiety, and fear) are controlled by braincentric neuromodulation systems (Guo, 2003). Previous studies on some cavefish species indicated that some brain genes may play an important role in the evolution of cave adapted behaviors. For example, a mutation in the MAO (monoamine oxidase) gene coding sequence (pro160leu) in the blind Mexican cavefish A. mexicanus brain may lead to high levels of serotonin and low monoamine oxidase activity, resulting in complete loss of aggressive behavior, allowing evolution of new food-seeking behaviors. The author believes that this is conducive to energy acquisition by A. mexicanus (Elipot et al., 2013(Elipot et al., , 2014. This may illustrate that the link between phenotypic, physiological and behavioral adaptations and their underlying genetic mechanisms in cave-dwelling fishes can be further understood by investigating the transcriptional regulation patterns of genes in cavefish brains. In this study, to minimize the influence of nature environment variables as much as possible, we selected two sympatric Sinocyclocheilus species, which represent an extremely rare phenomenon in the Sinocyclocheilus genus, namely, different types of fish (cave and surface) living in the same water body but at different depths. Through comparative transcriptomic analysis of brain transcriptome data which acquired through next-generation sequencing, followed by data analyses, we FIGURE 3 | The significant enrichment of DEGs KEGG pathway. The left of the X-axis shows the enrichment ratio, and the right bar shows the P value (p value is represented by color). *P < 0.05; **P < 0.01. FIGURE 4 | The number of differentially expressed genes of surface-cave fish groups. X-axis shows the differences in expression given as the log values, and the Y -axis shows significant differences in expression as FDR values in the Volcano plot; up-regulated genes and down-regulated genes are indicated by red dots and blue dots, respectively. identified tens of DEGs reliable candidate DEGs that may be related to Sinocyclocheilus cavefish adaptability from among 3,289 total DEGs.

DEGs That Involved in the Development of Sinocyclocheilus Cavefish Brain and Eyes
Compared to Sinocyclocheilus surface fish species, absent or great vestigial of eyes is one of the most remarkable troglodyte features in cavefish (Yoshizawa et al., 2012). Histological investigations of the eye tissue from Sinocyclocheilus cavefish have revealed that the cone cells normally distributed in the center of the retina are deficient, resulting in ineffectively distinguishment of colors and light intensity (Li and Tao, 2002). Another previously study of Sinocyclocheilus eyes have showed that gene expression patterns were substantially altered between cave and surface fish species, and some genes were identified that are responsible for eye degeneration (Meng et al., 2013). Interestingly, in this study, we also identified some genes (cryz[crystallin zeta], The left of the X-axis shows the enrichment ratio, and the right bar shows the P value (p value is represented by color). *P < 0.05; ***P < 0.001.  (Maeder et al., 2019). LCA symptoms are very similar to those of eye degeneration in Sinocyclocheilus cavefish. Unexpectedly, we found in this study that two LCA-related retinal developmental genes (CEP290 and nmnat1) were down regulated in the Sinocyclocheilus cavefish brain. CEP290 protein is encoded by 54 exons, and localizes to the photoreceptor-connecting cilium. Loss-of-function mutations in CEP290 can cause severe cone dystrophy and blindness (Cideciyan et al., 2007). nmnat1 encodes an essential rate-limiting enzyme that catalyzes the ligation of nicotinamide mononucleotide to ATP to form NAD in the NAD biosynthesis pathway, a highly energy-consuming process (Siemiatkowska et al., 2014). Earlier studies found that loss of nmnat1 in Drosophila melanogaster and mice leads to retinal degradation and photoreceptor cell degeneration (Zhai et al., 2006;Greenwald et al., 2016).
We believe that the DEGs related to eye development may reflect functions within the brain that regulate formation of the visual system of Sinocyclocheilus cavefish. We speculate that the differential expression of these genes may be one reasons for developmental defects in the eyes of Sinocyclocheilus cavefish eyes.
Among the Sinocyclocheilus cavefish species is a group of species with special head morphological characteristics. Namely, the head protrudes forward to form a spike-shaped angular convex structure. The function of this horn remains unknown. Li et al. speculated that the horn structure in cavefish species might be a protective organ that minimizes damage to the head when hitting underwater rocks (Li and Tao, 2002). In this study, we found two genes (coasy[Coenzyme A synthase] and pqbp1), which may regulate to head shape, that were significantly down regulation in the cavefish brain. CoA synthase, which is mainly expressed in the brain, is an important players in cellular metabolism and signal transduction. Loss of function variants of coasy have been shown to lead to head deformities (microcephaly) (Nemazanyy et al., 2006;van Dijk et al., 2018). Moreover, Khatri et al. (2016) found that down-regulation of coasy could lead to abnormal development of the head in zebrafish (Danio rerio), including brain protrusion in the knockout. Polyglutamine tract-binding protein 1 (PQBP1) has been directly linked to intellectual disability. Many studies have proven that mutations in PQBP1 (duplications or deletions of AG dinucleotides in the fourth coding exon) may be a primary cause of developmental deformities of the eyes (microphthalmia) and head (microcephaly) (Kleefstra et al., 2004;Martínez-Garay et al., 2007). Yang et al., 2020 analyzed the bone phenotype of Pqbp1-knockout mice (cKO) and found that Pqbp1-cKO mice exhibited brachycephaly and significant reduction in bone formation and chondrocytes. Furthermore, expression of osteoblast-and chondrocyte-related gene was downregulated in the Pqbp1-cKO mice. PQBP1, an intellectual disability causative gene, affects bone development and growth (Yang et al., 2020).
The DEGs involved in Sinocyclocheilus cavefish eye and head formation identified in this study may provide an appropriate genetic basis for future studies of the formation of troglomorphic traits and adaptability characteristics of troglobites. However, further validation experiments, such as gene knockout or gene overexpression, are needed to validate the function of these genes.

Downregulation of Insulin-Related Genes in the Sinocyclocheilus Cavefish Brain
Living in perpetual darkness in cave environments that cannot support photosynthesis have caused scarcity of food resources to become the main pressure faced by cavefish (i.e., adaptation to long periods of nutrient deprivation in underground water bodies). In addition to minimizing unnecessary consumption of energy (Gross and Wilkens, 2013), appropriate behavioral and metabolic adaptations to increase the probability of energy intake and improve maintenance of energy stores, such as blood sugar homeostasis and reduction of metabolic rate, are also important for survival.
Previous studies of the cave adaptability of fish have been found that, compared to surface fish, cavefish species Mexican tetra (A. mexicanus) exhibit dysregulated blood glucose homeostasis, and are insulin-resistant due to a mutation in the insulin receptor. The author believes that keeping blood sugar at high levels high for a long periods may be advantageous in caves as a strategy to cope with the periodic food shortages (Riddle et al., 2018).
Insulin is an essential polypeptide hormone that regulates of blood glucose, and fatty acid, and amino acid metabolism. It is synthesized exclusively in pancreatic islet B cells. Expression of the insulin gene is regulated by several islet-enriched transcription factors (Matsuoka et al., 2004). As part of a cascade of regulatory genes, the Maf transcription factor family has been shown to play an important role in insulin regulation.  7 | The annotation of the hub genes with top 30 connectivity in "blue" and "turquoise" module with differential expression identified by WGCNA.

Unigene id
Gene log2| FC| GO KEGG factor complex that activates insulin gene expression in pancreatic β-cells (Kataoka and Kamber, 2004), Artner et al. (2007) found that MafB loss of function mutations in mice lead to decreased insulin expression, and blockade of β-cells maturation. Similarly, Wataru et al. found that the number of pancreatic islet β-cells was significantly reduced upon knockout of the MafB gene in mice (Nishimura et al., 2006). MafA is a transcription factor that binds to the insulin gene promoter and is expressed exclusively in pancreatic β-cells to regulate insulin transcription in response to serum glucose levels.

Regulation of insulin
The current study demonstrated that the transcription factor MafA is a crucial regulator of insulin secretion. Moreover, the MafA-deficient mice are currently used as a new model for studying the overt diabetes mellitus (Matsuoka et al., 2004;Zhang et al., 2005).
In this study, we found that most of the Maf transcription factors family genes mentioned above, including MafA, MafB, and Maf K all trend toward down-regulation in Sinocyclocheilus cavefish species relative to surface fish species. Moreover, pax6 [paired box 6], the upstream regulatory gene of MafB and MafA, is synchronously down-regulated, though its log 2|Foldchange| < 1. Pax6 acts upstream of MafB, and is fundamental for MafB to activate insulin expression. Pax6 also regulates PDX-1 [pancreatic and duodenal homeobox 1], and MafA expression is required for β-cell maturation (Nishimura et al., 2008). We also found that, in addition to the Maf transcription factor family, other insulin secretion regulating genes, such as BRSK[BR-serine/threonine kinase-like protein], CDK16, and UBR1[ubiquitin protein ligase E3 component n-recognin 1] are differentially expressed in the cavefish brain. Brain-selective kinase (BRSK) is a member of the AMP-activated protein kinase (AMPK) subfamily of serine/threonine kinases. BRSK plays an important role in regulating glucose-stimulated insulin secretion (GSIS) in response to blood glucose by phosphorylating CDK16 (serine/threonine kinase) (Chen et al., 2012). UBR1 encodes one of at least four functionally overlapping E3 ubiquitin ligases of the N-end rule pathway. Previous studies have found that low expression of UBR1 perturbs pancreas cells and eventually leads to pancreatic insufficiency (Hülskamp et al., 2005). This may indicate that regulation of insulin expression to balance metabolic stress and energy intake is important in cavefish populations, but various mechanisms might be involved in different species. However, due to the limited fish sample sizes, we cannot conduct more additional experiments. Further functional verification experiments are needed to test the hypothesis that the Maf transcription factors family and other candidate genes we mentioned above play an important role in insulin regulation in cavefish is required.
Co-expressed Gene Network Revealed the Important Differences Between Two Sinocyclocheilus Fish Species In this study, we aimed to find "modules" which correlated significantly with different types of Sinocyclocheilus fish, by using transcriptomic data through WGCNA. Genes with certain functions in these modules may be related to the division of two ecological types of Sinocyclocheilus fish species. We focused on the two most highly correlated modules, "turquoise" and "blue, " and found three co-expressing network modules that comprised the top 90 hub genes (Figure 8). Through function and pathway annotation, we found that these hub genes with differential expression between cave and surface fish are mainly related to regulation of insulin, energy, and metabolism, phenotype and behavior, and nervous system development ( Table 7). Interestingly, most of the genes related to insulin regulation display varying degrees of differential expression in the brain of Sinocyclocheilus cavefish. This may further confirm that there is a difference in insulin regulatory mechanism between the two types of Sinocyclocheilus fish. Our previous research showed that two types of Sinocyclocheilus fish species have diverse dietary preferences and predation patterns (Chen et al., 2019). We believe that dietary preferences may cause differences in energy intake, ultimately affecting blood glucose regulation. Furthermore, these results may confirm the previous hypothesis that starvation resistance and metabolic regulation of cavefish would be present in multiple independently evolved cavefish populations (Bilandžija et al., 2020). However, more in-depth studies involving quantitative detection of insulin and blood glucose in the Sinocyclocheilus cave and surface fish species is still needed to verify the relationship between these DEGs and cave adaptability.
In this study, we conducted brain comparative transcriptomics analysis of two Sinocyclocheilus fish species, and found tens of DEGs that may be related to cave adaptability, with functions including insulin secretion regulation (MafA, MafB, MafK, BRSK, and CDK16) and troglomorphic trait formation (CEP290, nmnat1, COASY, and pqbp1). This study adds to our understanding of cave adaptation mechanisms. However, the brain transcriptome alone cannot fully represent and comprehensively explain the molecular mechanisms of cavefish special phenotypes. Accordingly, it is necessary to perform additional transcriptome sequencing from muscle, eyes, and liver in future work to provide a more complete explanation of the evolution of troglomorphism.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI BioProject (accession: PRJNA655363).

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Ethics and Welfare Committee at Yunnan University.

AUTHOR CONTRIBUTIONS
SC and HX contributed to conceptualization, project administration, and funding acquisition. CL, HC, and YZ carried out investigation and performed formal analysis. YZ and CL contributed to methodology and resources. HC, SC, and HX wrote the manuscript. HC and SC reviewed and edited the manuscript. All authors read and approved the final version of the manuscript.

FUNDING
This study was partially supported by the National Natural Science Foundation of China (31560111) and the Top Young Talents Program of Ten-Thousand Plan of Yunnan Province (YNWR-QNBJ-2018-024).