Integrated mRNA and microRNA transcriptomic analysis reveals the common and individual responses of zebrash embryos exposed to four types of graphene quantum dots (GQDs)

Background: Graphene quantum dots (GQDs) have great potential for bioimaging, biosensor, drug carrier, theranostics, and are recently reported as therapeutic agents to treat amylosis and inammation. Most types of GQDs have proven low toxicity in previous studies, but data about the transcriptomic responses of in vivo systems exposed to various GQDs remains largely unknown. Results: We examined the mRNA and microRNA (miRNA) expression changes of zebrash embryos exposed to four types of GQDs at a safe concentration (100 µg/mL), respectively as raw graphene quantum dots (R-GQDs), graphene oxide quantum dots (GOQDs), carboxyl GQDs (C-GQDs), and aminated GQDs (A-GQDs). The four GQDs elicited the number of differentially expressed genes (DEGs) in a decreasing order of A-GQDs, GOQDs, C-GQDs and R-GQDs to act on protein folding, potassium (K + ) and calcium (Ca 2+ ) channels, and spliceosome to varying degrees. Among the four GQDs, A-GQDs caused more genotoxic effects associated with lipid and hormone metabolism, MAPK signaling pathway, complement system, and ferropotosis. miRNA-seq data revealed that GOQDs aroused more differentially expressed miRNAs (DEMs), far exceed the total number of DEMs induced by the other three GQDs. dre-miR-735-5p and its potentially interactive gene-myogenin (myog) were identied as the only negatively-correlated miRNA-target gene pair shared by the four GQDs treatments. Conclusion: Taken together, this study provided substantial data underlying the common and specic transcriptomic responses of in vivo systems exposed to various types of GQDs, and also indicated the potential medicinal values of GQDs for ion antagonists and spliceosome-targeted therapies.


Background
Graphene quantum dots (GQDs) are newly-emerged graphene-derived nanoparticles (NPs) with lateral dimensions typically below 10 nm in sizes, and show better performance in electronic, optical, spin, and photoelectric properties induced by the quantum con nement effect and edge effect [1]. GQDs hold promise in biomedical applications for their high biocompatibilities [2], tunable photoluminescence (PL) [3], light-dependent anti-(pro-) oxidant activities [4], and intrinsic peroxidase-like catalytic activities [5,6], which have been exploited as sensitive uorophores for real-time bioimaging and anticancer photodynamic therapy in in vivo systems [7][8][9], and also independently used for anti-in ammation via their strong free radical scavenging activities [10][11][12]. Interestingly, two recent studies found that GQDs could speci cally disaggregate amyloid proteins (i.e. α-synuclein brils, islet amyloid polypeptide (IAPP)) and reduce the neurological toxicity in Parkinson's disease as well as improve IAAP-related symptoms [13,14]. Although most studies claimed that GQDs are well biocompatible and induce low toxicity in organisms, there remain a few inconsistent ndings regarding the potential toxicity of GQDs that may result from nanofabrication methods and the varying surface characters, including nanosizes, surface charges, functional groups, element doping, and impurities [15,16]. It is now generally believed that hydroxylated GQDs (OH-GQDs) are hazardous for its inducibility of oxidative stress even at relatively low exposure concentrations, leading DNA damage, G 0 /G 1 cell cycle arrest, the disruption of microtubule structure, enhanced intestinal penetrability, and activation of multiple signaling pathways [17][18][19]. At present, nitrogen-doped GQDs (N-GQDs) and aminated GQDs (A-GQDs) have been investigated much for their toxicity, although some previous studies asserted that they were low toxic to cell survival and animal development [18, [20][21][22]. Most toxic data about the two GQDs were evidenced from their effects on cellular micromorphology, nucleotide structures, and gene expression pro ling. Both were found to induce DNA methylation and transcriptomic responses in relation to the canonical MAPK signaling pathway [23][24][25][26]. In addition, it was reported that N-GQDs could deform the structure of lipid droplets, disrupt the redox-sensitive system, and cause ferroptosis [27][28][29], while A-GQDs could induce DNA cleavage via Hbonding and π-π stacking [30]. In general, the low toxicity of most types of GQDs was concluded from the conventional toxicity testing, involving the assessment of cell growth and survival, the developmental indices in animal models, and the canonical blood/molecular markers [22,[31][32][33], and there remains incomplete whether all GQDs share the common toxic pathways and to what extent various types of GQDs perturb speci c responses of genes, proteins, metabolites, as well as the behaviors. Now there have been reported several studies as mentioned above that employed toxicomics approaches to unveil the potential genotoxicity caused by GQDs [24][25][26], but the relevant toxicity assessment of GQDs using highthroughput screening methods based on in vivo systems are seldom [24].
Zebra sh as an ideal animal model has been extensively used in the biosafety evaluation of a wide range of engineered nanomaterials (ENMs) [34]. This animal possesses a suit of conserved evolutionary system that is comparable to that in vertebrates and has developed a typical optical-transparency body at early stage that is easily observed when undergoing phenotypic variations caused by genetic changes, environmental harmful factors, and among others [35]. Compared to the mouse model, zebra sh exhibits unique characteristic in its high-fecundity that provides considerable number of embryos for highthroughput drug screening [36]. In addition, the sh is low-cost reared and easily genetically manipulated in sizes. Zebra sh has been widely applied for the toxicity assessment of a series of carbon-based nanomaterials, including carbon nanotubes (CNTs) [37], carbon fullerene (mainly C60) [38], graphene (G) [39], graphene oxide (GO) [40], carbon quantum dots (CQDs) [41], and N-GQDs [28,42]. In the study, we employed high-throughput sequencing technology to explore the potential side effects of four types of GQDs on in vivo transcriptional pro ling as opposed to gene selection based on known markers in canonical toxicity assays. The study aimed to reveal the common gene clusters and associated signaling pathways induced by the four GQDs in in vivo systems, and also pointed out the speci c pathway and the miRNA-mRNA interaction triggered by different GQDs. This study would provide valuable insights about the potential genotoxicity of various GQDs on in vivo systems, and also an alert for the selection of GQDs especially for pharmacotherapies.

Zebra sh Husbandry And Exposure Assays
From the larval stage, a batch of wild-type zebra sh (AB) were raised according to the standard breeding protocol [43]. Zebra sh at about 4 months old were coupled for spawning and the newly-laid eggs were collected at four post hour fertilization. Fifty embryos were slightly pipetted into every well of 6-well cell plates lled with 3 mL various GQDs solutions (100 µg/mL), which was newly prepared every day. Each treatment was performed for at least four biological replicates, and control group was only treated with water. The survival rates of zebra sh embryos exposed to different GQDs were recorded at 2 day post fertilization (dpf) until to the end of the assay. Hatching rate was estimated at 2 dpf and 3 dpf, and heart rate (beats/min) was observed at 3 dpf, 4dpf, and 5dpf using an inverted microscope (Nikon SMZ18, Japan), each treatment was carried out with at least three independent biological replicates. After exposure for 7d, zebra sh embryos were collected from each well and kept at -80℃ refrigerator.

MRNA And MIRNA Extraction And Sequencing
Total RNAs were extracted from zebra sh embryos using TRIzol® Reagent according to the manufacturer's instruction (Invitrogen). Genomic DNA debris was degraded using DNase I RNase-free (Takara, Japan). RNA quality was evaluated on 1% agarose gels and its concentration was measured by means of Nanodrop-2000 (Thermo Fisher Scienti c, USA). Finally, RNA integrity was monitored using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and only high-quality RNA samples (OD260/280 = 1.8 ~ 2.2, OD260/230 ≥ 2.0, RIN ≥ 7, 28S:18S ≥ 1.0, quantity > 5µg) were applied for constructing sequencing libraries. (1) mRNA sequencing: RNA-seq transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA) using 1µg of total RNA. Shortly, messenger RNA was isolated according to polyA selection method by oligo (dT) beads and then fragmented by fragmentation buffer rstly. Secondly, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end-repair, phosphorylation and 'A' base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR ampli ed using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quanti ed by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten/NovaSeq 6000 sequencer (2 × 150bp read length). (2) miRNA sequencing: RNAseq transcriptome library was prepared following TruSeq™ RNA sample preparation Kit from Illumina (San Diego, CA) using 1µg of total RNA. Shortly, messenger RNA was isolated according to polyA selection method by oligo(dT) beads and then fragmented by fragmentation buffer rstly. Secondly double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end-repair, phosphorylation and 'A' base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR ampli ed using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quanti ed by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten/NovaSeq 6000 sequencer (2 × 150bp read length).

Bioinformatics Analysis
Principal component analysis (PCA) was performed for mRNAs and miRNAs data with OmicShare tools, a free online platform for data analysis (http://www.omicshare.com/tools). The fragments per kilobase million (FPKM) and transcripts per million (TPM) values were employed to quantify mRNA and miRNA expression levels, respectively. The differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) were screened out using DESeq2. Gene ontology (GO) and KEGG pathways of DEGs were analyzed with the online DAVID Bioinformatics Resources 6.8 and KOBAS 3.0 software, with a Bonferroni-corrected P-value ≤ 0.05 compared with the whole-transcriptome background. Heatmap of DEGs and DEMs among different treatment was drawn using ggplot2. Protein-protein interaction (PPI) network was deduced by means of online String tool and Cytoscape (Version 3.8.0). Through a BLAST search of the miRbase (version 21.0), the perfectly matched sequences were used to count and analyze the known miRNA expression pro le. In-house scripts were used to obtain the identi ed miRNA base bias on the rst position with certain length and on each position of all identi ed miRNA. The miRNA-target interactions were predicted by OmniSearch analysis that provides access to multiple target prediction dataset, including TargetSCan, miRanda, miRDB, and RNAhybrid. An index of corrected correlation < − 0.5 (Spearman's correlation coe cients) was used to evaluate the negative correlations between known miRNA and target genes.  Table S1. The relative mRNA and miRNA levels were respectively normalized to the internal control egfα and U6 snRNA and calculated using the 2 −ΔΔCt method. Data were shown as Mean ± SD from four independent biological replicates.

Results
Surface characters and developmental toxicity assessment Transmission Electron Microscopy (TEM) was used to characterize the morphology and size distribution of four kinds of GQDs. As shown in Fig. 1a, the TEM images showed that the four kinds of GQDs gave a narrow size distribution with well monodisperse, which indicated their high quality. Clear lattice fringes exhibited their graphitic structure with lattice spacing (d = 0.21nm, 0.24 nm) corresponding to (100) diffraction facets of graphite [44][45][46]. Data of XPS full spectrum scan showed that GOQDs were among the most complex in the atomic ratio, including C1s (44.21%), O1s (34.59%), Na1s (6.79%), Si2p (8%), Cl2p (2.44%), S2p (2.28%), and Sc2p (6.11%), whilst C-GQDs contained only C1s (58.13%) and O1s (41.87%) ( Fig. 1b and Table S2). We performed the curve tting analysis based on the individual elemental scanning, mainly C1s and O1s (Fig. S1a). C1s spectra of R-GQDs and A-GQDs were deconvoluted into ve peaks of C-C/C = C (~ 284. eV), and NH n + (400.8 eV) (Fig.S1b), which were consistent with the previous report [30]. Next, it was found from the zeta potential analysis that the four GQDs were all negatively charged in three mediums (Fig.S1c), including ddH 2 O, 1 × PBS, and 1 × DMEM. Among the four NPs, C-GQDs induced the weakest zeta potential changes in the three mediums, ranging from − 1~ -5 eV. The developmental characters of zebra sh embryos including survival rate, hatching rate, and hearbeat rate were not signi cantly changed by different GQDs treatment (100 µg/mL) when compared to control group ( Fig. S2: a-c). Our exposure assay manifested the low toxicity of the four GQDs to zebra sh embryonic development.

Whole Transcriptomic Pro ling
Principle component analysis (PCA) revealed that the mRNA expression data extracted from GOQDs and C-GQDs treatment clustered in close distances, as well as the known miRNA sequencing data among the clusters of control, R-GQDs, and C-GQDs treated groups. A-GQDs treatment was clearly separated from the others regardless of mRNA or miRNA data ( Fig. 2a and 2b), and the miRNA data from GOQDs treatment also clustered independently. With the same corrected p value (FDR) < 0.05 and two cutoff settings (FC > 2 or 1.5), two ensembles of DEGs were screened out between control and NPs treated groups. A-GQDs elicited the largest number of DEGs (FC > 2, 1,878; FC > 1.5, 4,796), far more than the total number of DEGs induced in other three GQDs exposed groups (Fig. 2c, Table S3). However, from the miRNA expression levels, GOQDs were found to induce more DEMs (FC > 2, 32; FC > 1.5, 65) than other GQDs, and R-GQDs and C-GQDs exerted negligible effects on miRNAs expression levels (Fig. 2c, Table  S4). To excavate the potential genes commonly induced by different GQDs, the Venn analysis based on the DEGs and DEMs data from respective GQDs treatment was performed with a FC cutoff > 1.5 and a signi cant FDR value of < 0.05, and 75 DEGs were identi ed to be shared by the four GQDs treatments (Fig. 2d, Table S5), and no overlapping miRNAs are found among various GQDs treatments.

Data Pro ling Of DEGs
Conventionally, we performed the functional enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for DEGs identi ed from different GQDs treatment. With the setting of the threshold as FDR < 0.05, more than 500 subcategories were signi cantly enriched into three divisions in A-GQDs treated group, encompassing biological process (BP, 329 terms), cellular component (CC, 86 terms), and molecular function (MF, 91 terms). Comparatively, the signi cantly-enriched GO items were less found in C-GQDs (66), GOQDs (48), and R-GQDs treated group (48) ( Table S6). In the top ten GO terms from A-GQDs and GOQDs groups, both share two identical GO terms -'cellular response to unfolded protein and response to unfolded protein' (Fig. 3a), from which DEGs related to hsp genes are highly upregulated. Compared to other GQDs treatments (Fig. 3b), A-GQDs treatment induced more KEGG pathways in relation to lipid and hormone metabolism (Fig. S3a), MAPK signaling pathway (Fig. S3b), complement system (Fig. S3c), ferropotosis (Fig. S3d), cytochrome P450 system (Fig. S3e), all of which directly or indirectly interacted with each other and nally triggered the (pro-) apoptotic signals (Fig. S4).
The up-regulation of hsp genes usually signals the organisms suffering stressful environmental conditions and their translated proteins (heat shock proteins, HSPs) act as intra-cellular chaperones to stabilize novel proteins by guiding their proper folding and conformation [47], we herein made a generalization for all hsp genes and relevant genes across the mRNA expression pro ling in all treatments. Heatmap displayed that A-GQDs induce the most number of hsp and associated genes, and R-GQDs elicited the minimal transcriptomic responses of HSP family members (Fig. 4a). To validate the RNA-seq data, ten hsp genes were selected from the Heatmap for qPCR ampli cation. Results con rmed that most of the qPCR data were in accordance with the mRNA-seq results, and hspa4a, hsp70.2, hsp70.3, and hsp90aa1.2 were signi cantly upregulated in the four GQDs treatments when compared to controls (Fig. 4b).
To better explore DEGs and associated signaling pathways commonly aroused by the four GQDs, we rstly performed a Venn analysis based on the signi cantly enriched GO terms from different GQDs treatment, and two GO terms -'ion membrane transport and ion transport' are found to be shared by the four GQDs treatments (Fig. 5a). The further Venn analysis revealed that 6 and 10 DEGs were shared by the four GQDs in each GO term ( Fig. 5b and Fig. 5c), and nally 8 genes were determined from these DEGs (Fig. 5d). Except the upregulated expression of apoa4b.2 in the four GQDs treatments, other genes related to potassium (K + ) and calcium (Ca 2+ ) channels were signi cantly suppressed by the four GQDs. We further made a generalization about the overlapping 75 DEGs among the four GQDs treatments, which exhibited a congruous expression pro ling within three repeating samples between control and treatment groups (Fig. 5e). GO enrichment analysis showed that the overlapping genes were signi cantly enriched into protein N-acetylglucosaminyltransferase and RNA processing, and splicesome was identi ed as the only signi cantly enriched KEGG pathway (Table S7). Next, we performed the proteinprotein interactions (PPI) analysis based on these overlapping genes and showed that serine/argininerich splicing factor 11 (srfs11) as the center mediator coordinating with other splicing factor to regulate downstream genes, including nktr, hspbp1, rbm5, gpatch8, and arglu1a (Fig. 5f).

Data Pro ling Of MIRNAs
In the study, the target genes of DEMs in different GQDs treatment were screened out using multiple prediction databases with reference to respective DEGs database, from which a total of 148 and 153 coexpressed genes for 52 and 13 miRNAs were respectively identi ed from GOQDs and A-GQDs treatments, and was few miRNA-mRNA pairs found in the other two treatments. Owing to very few GO and KEGG pathways signi cantly enriched from the target DEGs (FDR < 0.05), we focused on the data mining of the miRNA-mRNA pairs with negative correlations from the four GQDs treatments. 86 and 42 anti-correlated miRNA-mRNA pairs were identi ed from GOQDs and A-GQDs treated groups, respectively extracted from the combination of 43 miRNAs: 65 genes in GOQDs and 11 miRNAs: 36 genes in A-GQDs treatment, while only one negative miRNA-mRNA pair was identi ed from the other two treatments (Table S8). Together with the above generalization of the common signaling pathways induced by the four GQDs, several pairs from GOQDs likely involved in the ion membrane transport, ion transport, and the spliceosome pathway. From the GOQDs treatment, three genes (slc25a23a, kcnn1a, trpm1b) implicating in the 'ion membrane transport and ion transport' were predicted to be repressed by four miRNAs (dre-miR-455-2-5p, dre-miR-140-5p, 10_1985, 24_18026), while three genes (srsf2a, srsf5a, srsf5a) involving in the spliceosome pathway were probably negatively regulated by dre-miR-15a-5p and dre-miR-16a. In addition, three genes (socs1a, il6st, map2k6) involving in the MAPK signaling pathways were identi ed to be repressed by dre-miR-7132-3p, dre-miR-19a-5p, and dre-miR-216a in GOQDs treatment. dre-miR-29b3-3p from A-GQDs treatment showed the strongest multitargeting activities and likely involved in the downregulated the mRNA levels of nine genes, in which two pairs (dre-miR-29b3-3p: slc38a8b and dre-miR-29b3-3p: slc44a2) may collaboratively act on the splicesome pathway.
From all negatively correlated miRNA-mRNA pairs, we found that the pair of dre-miR-735-5p and its negatively regulated gene -myogenin (myog) was shared by A-GQDs, GOQDs, and C-GQDs treated groups, and dre-miR-735-5p was not signi cantly up-regulated by R-GQDs according to the originally-used screening threshold (FDR < 0.5, FC > 1.5). However, it was found that dre-miR-735-5p was induced more than three-fold change (3.165) in R-GQDs treated group when the cutoff as set as p value < 0.01 and FC > 1.5. Then, we employed qPCR to investigate the expression changes of dre-miR-735-5p and its target gene-myog among the four GQDs treatments, and results showed that the expression levels of dre-miR-735-5p in all GQDs treatments were signi cantly activated when compared to control group, and myog exhibited the similar expression pro ling to that in RNA-seq data (Fig. 6), and then suggesting that the dre-miR-735-5p-myog pair could be commonly aroused by the four GQDs.

Discussion
GQDs are typical of graphene structure with honeycomb lattices and inherit many physicochemical properties from graphene. GQDs develop strong photoluminescence (PL) and photostability with highly tunable band gaps, which have been exploited for sensitive uorophores to monitor living cells and drug carriers to assist target therapy in in vivo systems. More strikingly, several recent studies showed that GQDs can be independently used as pharmaceutical agents to treat in ammations and amylosis [13,14], thereby being paid much attention for its potential medicinal value. However, GQDs are diversely characterized by the fabrication method and raw material, thus leading to the potential nano-bio effects that may be unsafe [48]. Therefore, the systemic biosafety assessment of GQDs would be demand before its approval for clinical trials. In the previous study, we had examined the transcriptomic responses (mRNAs) of zebra sh embryos exposed to N-GQDs, which were found to trigger pathways mainly associated with acute in ammation and detoxifying process [24]. In the study, we extended the number of GQDs for toxicity evaluation, and the four types of GQDs including R-GQDs, GOQDs, C-GQDs, and A-GQDs were the rst time to be investigated for their potential genotoxic effects based on in vivo systems. Our exposure assay exhibited the low toxicity of the four GQDs on the developmental characters of zebra sh embryos, some of which was also observed in mice treated with the same four GQDs [22], and then suggesting that the adverse effects induced by GQDs would not be easily discriminated upon the conventional toxicity assays such as the assessment of cell growth and survival, the developmental indices in animal models, and the canonical blood/molecular markers. Thus, the employment of highthroughput transcriptomics would be an optimal choice to unravel the potential genotoxic effects induced by various types of GQDs. Our study showed that the four GQDs can act on unfolded protein responses and endoplasmic reticulum via activating the HSPs family members and associated pathways to varying degrees, then suggesting that GQDs may cause dysfunction of certain proteins, but the process of which would be counteracted by a group of HSP family members mobilized by GQDs. Previous studies showed that GQDs can break down amyloid oligomerization and are promised as potential small molecules for treating amylosis [13,14], and then whether GQDs speci cally target amyloid proteins or not remains to be explored in the future.
Human body is driven by a complex electric circuit comprising nely tuned and intricately synchronized charge transfer systems, which are the principal drive for the ows of energy to sustain the vital activities, such as the Krebs cycle (citric acid or tricarboxylic acid cycle) in aerobic organisms, being a singledirection cyclic pathway transferring electrons from citrate to NAD + , FAD and GDP to produce NADH, FADH 2 , and GTP [49]. Importantly, the electron transferring serves as bidirectional messengers in the redox homeostasis, which maintains the normal physiological steady state in living organisms [50]. GQDs show metal-like attributes and have been reported to possess anti-and pro-oxidant activities [4], and our functional enrichment analysis revealed that the four GQDs can disturb K + and Ca 2+ homeostasis by co-suppressing K + and Ca 2+ channel associated gene clusters (Fig. 5C), then suggesting that GQDs may involve in the metal homeostasis via energy and charge transfer. The inhibitory effects of GQDs on K + and Ca 2+ channel were also observed in the mouse lung and the primary cortical neurons of rats treated with graphene oxide nanosheets [51,52]. Blocking K + and Ca 2+ channel can jointly or independently delay the action potential repolarization and intensify the blood vessel dilatation, and sometimes trigger irregular physiological responses in muscle and cardiovascular system [53,54]. However, this impeding effect has been often reported to play active roles in treating arrhythmia [55], hypoglycemia [56], and hyperglycemia [57], thus meaning GQDs could be potential ion channel antagonists to treat relevant diseases. Therefore, the inhibitory effects of GQDs on K + and Ca 2+ channel should be adequately evaluated by combing with other physiological responses. Our nding about the transcriptional regulation of spliceosome triggered by GQDs was also reported in A549 cell line and zebra sh treated with GO [58,59], implying that the pathway could be commonly aroused by GQDs and other graphene derivatives. Splicing factors such as splicing factors 3b (SF3B) have been reported as a core target for anticancer treatment [60], and in this study ve splicing factors (srsf2a, sf3b1, srf4, srsf5a, and zgc: 163098) from the spliceosome pathway were all signi cantly downregulated by the four GQDs, indicating the potential of GQDs as spliceosome inhibitors or adjuvants for targeted cancer therapy.
In this study, the number of miRNAs triggered by GOQDs was much more than that of other three GQDs.
However, the miRNA clusters and associated miRNAs-mRNAs regulatory network induced by GOQDs were totally different from that aroused by GO in in vivo system [61,62], suggesting that the regulatory effects of GOQDs on miRNAs have been developed independently from that of G and GO. Still, a few miRNA-mRNA pairs identi ed from GOQDs treatment may involve in MAPK signaling pathway (Table S7). Although the miRNA responses were weak relative to the mRNA responses induced by GQDs, a negatively regulated miRNA-mRNA (dre-miR-735-5p: myog) was found to be aroused by the four GQDs. Myogenin (myog) is reported as the key regulator for the proper differentiation of most myogenic precursor cells during the process of myogenesis [63], indicating that GQDs may preclude the development of functional skeletal muscle via the commonly activation of dre-miR-735-5p and its target gene-myog.
Among the four GQDs, A-GQDs showed more genotoxic effects on in vivo system compared to other GQDs, and most of the toxic pathways aroused by A-GQDs were also identi ed in N-GQDs treated zebra sh embryos [24], including the complement and coagulation system, MAPK pathway, cytochrome P450 system, lipid and hormone metabolism, etc. Interestingly, we found that A-GQDs and N-GQDs shared the similar N1s pro le in two peaks corresponding with amide bond and ammonium ions [42]. Together with the previous nding about the disruptive effects of A-GQDs on DNA chains [30], the high genotoxicity of A-GQDs may result from the enhancement of the polarization of A-GQDs caused by the lone pair electrons of the amino groups, which react with the phosphoric acid in nucleotides and disturb the structure of DNA bases. Therefore, A-GQDs should be functionally improved to lessen its potential genotoxic effects when used for the in vivo therapies.

Conclusion
In the present study, we con rmed that the four GQDs at the same exposure concentration (100 µg/mL) had no noticeably adverse effects on zebra sh embryonic development. Our microarray data was the rst time to reveal that various types of GQDs act on at least three gene clusters in relation to protein folding, ion channels, and RNA splicing to varying degrees. A-GQDs were among the most genotoxic NPs for its perturbation of multiple signaling pathways, most of which were also identi ed from that in zebra sh exposed to N-GQDs. Although the miRNAs expression changes aroused by the four GQDs were much weaker than mRNAs responses, a negatively correlated miRNA-target gene pair (dre-miR-735-5p: myog) was identi ed to be commonly activated by the four GQDs, showing a common miRNA-mRNA regulatory network occurring among all GQDs. Our mRNA-seq data also provided invaluable information for the design of targeted drug delivery based on GQDs biocompatibility.