Integrated mRNA and miRNA transcriptome analysis provides novel insights into the molecular mechanisms underlying goose pituitary development during the embryo-to-hatchling transition

It is well established that the endocrine system plays a pivotal role in preparing the avian embryos for the abrupt switch from chorioallantoic to pulmonary respiration during the critical embryo-to-hatchling transition. However, as the master gland of the endocrine system, there has been little research focusing on the molecular mechanisms controlling the development and function of the pituitary gland during the peri-hatch period in birds. In the present study, we aimed to determine the genome-wide mRNA and miRNA transcriptome profiles of the pituitary during the embryo-to-hatchling transition period from embryonic day 22 (E22) to post-hatching day 6 (P6) in the goose (Anser cygnoides). Of note, expression of Anser_cygnoides_newGene_32456 and LOC106031011 were significantly different among these 4 stages (i.e., E22, E26, P2, and P6). Meanwhile, the neuroactive ligand-receptor interaction pathway was significantly enriched by the DEGs commonly identified among three pairwise comparisons. At the miRNA transcriptome level, there were not commonly identified DE miRNAs among these 4 stages, while the 418 of their predicted target genes were mutually shared. Both the target genes of DE miRNAs in each comparison and these 418 shared target genes were significantly enriched in the ECM-receptor interaction and focal adhesion pathways. In the predicted miRNA-mRNA interaction networks of these 2 pathways, novel_miRNA_467, novel_miRNA_154, and novel_miRNA_340 were the hub miRNAs. In addition, multiple DE miRNAs also showed predicted target relationships with the DEGs associated with extracellular matrix (ECM) components. Among them, expression of novel_miR_120, tgu-miR-92-3p, and novel_miR_398 was significantly negatively correlated with that of LAMC3 (laminin subunit gamma3), suggesting that these miRNAs may regulate pituitary tissue remodeling and functional changes through targeting LAMC3 during development. These identified DE mRNAs and miRNAs as well as their predicted interaction networks involved in regulation of tissue remodeling and cellular functions were most likely to play critical roles in facilitating the embryo-to-hatchling transition. These results provide novel insights into the early developmental process of avian pituitary gland and will help better understand the underlying molecular mechanisms.


INTRODUCTION
Pituitary is known as the master gland of the endocrine system in vertebrates and plays essential roles in a wide spectrum of biological activities such as growth and development, immune and metabolism, and reproduction (Musumeci et al., 2015). Anatomically, the vertebrate pituitary is comprised of multiple cell types, including corticotrophs, somatotrophs, lactotrophs, gonadotrophs, thytrophs, melanocyte stimulating hormone (MSH) secretory cells, and folliculo-stallate (FS) cells (Le Tissier et al., 2012). The orderly differentiation and functional maturation of these cell components are undoubtedly a prerequisite for ensuring the versatile functions of the pituitary (Ooi et al., 2004). Of note, it has been previously reported that both the anatomical structures and developmental dynamics of the pituitary gland vary significantly between mammals and birds 1 (Scanes, 2015), especially during embryonic and early post-hatching life (Scanes et al., 2005). Since the embryo's responses to new environmental challenges may have short-and long-term impacts on the physiology and performance, it is therefore of particular interest to dissect the molecular mechanisms regulating the development and function of the avian pituitary during the peri-hatch period, eventually helping better understand how the embryo adapts itself to the abrupt environmental changes during the critical embryo-tohatchling transition.
In order to adapt well to the sudden environmental changes during the embryo-to-hatchling transition, both the morphology and functions of liver, thyroid, and ovary have changed simultaneously during this period in birds (McNabb, 2006;Cogburn et al., 2018;Hu et al., 2020). As the critical regulator of the development and functions of these tissues or organs, there is evidence that the pituitary histomorphology and gene expression profiles also change remarkably during early developmental stage. Moreover, species-specific differences in gene expression profiles during early pituitary development have been previously shown in domestic birds. In the chicken pituitary, the mRNA expression levels of growth hormone (GH) and prolactin (PRL) started to increase on embryonic day 16 (E16) and E18, respectively (Kansaku et al., 1994), while those of duck GH and IGF1 started to increase on E25 and E21, respectively (Hu et al., 2015). Although the chicken pituitary transcriptome changes in the posthatching development have been previously investigated (Pritchett et al., 2017;Ellestad et al., 2019), there are still lack of studies on the avian pituitary developmental dynamics during the critical embryo-to-hatchling transition, especially in waterfowls.
In recent years, high-throughput sequencing technologies have been demonstrated to be efficient and reliable to study the developmental dynamics of the avian muscle and ovarian tissues during the embryonic to early post-hatching stages (Liu et al., 2019). In the meantime, as a class of non-coding RNAs, miRNAs have also been shown to be involved in regulating the growth, development, and metabolism of multiple avian tissues (Fu et al., 2018;Li et al., 2019). Therefore, the objectives of this study were to: 1) explore the dynamic changes in the goose pituitary mRNA and miRNA transcriptomes during the peri-hatch period and 2) identify the key mRNAs, miRNAs as well as the miRNA-mRNA interaction networks involved in facilitating the embryo-tohatchling transition. These data will provide new insights into the molecular mechanisms underlying early pituitary development in birds.

Ethics Statement
All geese were obtained from the Waterfowl Breeding Experimental Farm of Sichuan Agricultural University. All experimental procedures that involved in animal manipulation were approved by the Institutional

Sample Collection and RNA Extraction
Four representative developmental stages during goose embryo-to-hatchling transition, including E22, E26, post-hatching day 2 (P2), and P6, were selected for sample collection. Six female geese with similar weights at each developmental stage were randomly selected from 150 geese hatched from the same batch and euthanized by carbon dioxide anesthesia and exsanguination by severing the carotid artery, and the pituitary gland from each individual was snap-frozen into liquid nitrogen and finally stored at À80℃ until RNA isolation. According to the manufacturer's protocol, the total RNA was extracted from all collected samples using the miRNeasy Mini Kit (Qiagen, Hilden, Germany). The extraction, purity, concentration, and integrity of RNA were determined by Nanodrop (Thermo Fisher Scientific, Wilmington, DE) and Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). The criteria used to select the RNA for following analysis were A260/A280≥1.8, A260/A230≥2.0 and RNA integrity number >7.0.

Library Construction and Sequencing
To satisfy the requirements for sequencing, 2 hypophyses were pooled for each replicate on E22, E26, P2, and P6 (Supplementary Table S4). Twelve mixed RNA samples were used for library construction. Total RNA was used to prepare a small RNA (sRNA) library according to instructions of NEBNext Multiplex sRNA Library Prep Set for Illumina (Illumina, San Diego, CA). The sRNA library was validated by using an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA) to check its size and purity. The NEBNext Poly (A) mRNA magnetic separation module (NEB, E7490, Ipswich, MA) was used to isolate mRNA. The cDNA library was constructed using the NEBNext Ultra RNA Library Prep Kit from Illumina (NEB, E7530, Ipswich, MA) and the NEB Multiplex Oligos from Illumina (NEB, E7500, Ipswich, MA) as per the manufacturer's instructions. All libraries were sequenced using Illumina HiSeq 2500 (Illumina, San Diego, CA). The original sequencing data for this study can be found in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) at NCBI, with the BioProject ID: PRJNA649062.

Bioinformatic Analysis
Raw data of RNA-seq were processed using the NGS QC Toolkit (version 2.3.3) (Patel and Jain, 2012). Reads that contained poly-N and low-quality reads were removed to obtain clean reads. Clean reads were mapped to the goose reference genome (https://www.ncbi.nlm. nih.gov/assembly/GCF_000971095.1/) using the HISAT2 software (version 2.1.0) (Daehwan et al., 2015). StringTie (version 1.3.3b) (Mihaela et al., 2015) was used to assemble transcripts after integrating all individual transcripts and genes.
Raw data in fastq format were firstly processed using in-house perl scripts. Clean data were obtained by removing reads that contain adapters, ploy-N, and lowquality reads. Reads were trimmed and cleaned by removing the sequences smaller than 18nt or longer than 30nt. Bowtie software (version 1.2.3) was used for sequence alignment of clean reads and Silva databases, GtRNAdb databases, Rfam databases, and Repbase databases. The remaining reads were used to detect known miRNA and novel miRNA predicted by comparing with Genome. The potential target genes of miRNAs were predicted using miRanda (Betel et al., 2008) and targetscan (Lewis et al., 2003) softwares.
Differential expression analysis among groups was performed using the R package EBseq (version 1.29.0) (Leng et al., 2013). The P values were adjusted using the Benjamini and Hochberg's approach in order to control the false discovery rate. The RNAs and miRNAs with an adjusted P-value (FDR) <0.05 and |Fold change| >1.5 analyzed by EBseq were regarded as DE mRNAs and miRNAs. Gene ontology enrichment analysis software tools (GOEAST) (Zeng and Wang, 2008) was used to analyze the Gene Ontology (GO) functions. KOBAS3.0 (Ai and Kong, 2018) was used to analyze the Kyoto Encyclopedia of Genes and Genomes (KEGG) functions.
Cytoscape software (version 3.7.2) (Shannon et al., 2003) was used to show the targeting relationship between DE miRNAs and their target genes. String online software (Szklarczyk et al., 2019) was used to display protein interaction networks between target genes.

Quantitative Real-Time PCR Analysis
Total RNA extracted from each pituitary was reverse transcribed into cDNA using a Revert Aid First Strand cDNA Synthesis Kit (Thermo, MA). Primer 5.0 was used to design the primers for quantitative real-time PCR (qRT-PCR) ( Table 1). The PCR reactions were performed on the CFX96TM Real-Time PCR Detection System (BioRad, Hercules, CA) using the SYBR Premix Ex TaqTM II (Takara Biotechnology Co., Ltd., Dalian, China). Each sample was repeated in 3 times. The relative expression levels of the target genes in the samples were normalized to the reference genes GAPDH and b-ACTIN using the comparative Cq method (DDCq) (Schmittgen and Livak, 2008).

Statistical Analysis
The expression levels of DEG were expressed as the mean § SEM. Pearson's correlation coefficient was calculated to analyze the correlations between the expression of LAMC3 and several miRNAs (novel_miR_120, tgu-miR-92-3p, and novel_miR_398). Duncan's Table 1. Primers used for qRT-PCR in this study.

Primer name
Sequence (5' -3') Product length (bp) multiple range test was used to compare the differences in the expression levels of LAMC3, novel_miR_120, tgu-miR-92-3p, and novel_miR_398 at E22, E26, P2, and P6. P-value <0.05 was considered statistically significant. All statistical analyses were carried out using the SPSS 27.0 software.

Overview of the mRNA Transcriptome of Goose Pituitary During the Peri-hatch Period
A total of 99.57 Gb clean bases were obtained from 12 samples through mRNA sequencing, and 75.69%8 0.84% of the sequencing reads were aligned to the goose reference genome (Supplementary Table S1). As shown in Supplementary Figures 1A and 1B, 3 biological replications at each sampling point were well clustered together. Moreover, three biological replications on E22 were closer to those on E26, so were those on P2 and on P6. There were 303, 880, and 118 DEGs identified between E22_vs._E26, E26_vs._P2, and P2_vs._P6, respectively ( Figure 1A). Meanwhile, 74, 33, and 11 DEGs were overlapped in E22_vs._E26 and E26_vs. _P2, E26_vs._P2 and P2_vs._P6, and E22_vs._E26 and P2_vs._P6, respectively ( Figure 1B). As shown in Figure 1B, there were 2 genes (i.e., Anser_cygnoides_-newGene_32456 and LOC106031011) whose expression in each sampling point was significantly different.

Functional Enrichment Analysis of DEGs Between Goose Pituitaries at Different Developmental Stages
The top 20 KEGG pathways enriched by DEGs in 3 pairwise comparisons were shown in Figures 2A-2C. Of them, the AGE-RAGE signaling in diabetic complications, apelin signaling, adrenergic signaling in cardiomyocytes, and ECM-receptor interaction pathways were commonly enriched by the DEGs in E22_vs. _E26 and E26_vs._P2. The focal adhesion and MAPK signaling pathways were commonly enriched by the DEGs in P2_vs._P6 and E26_vs._P2. However, the calcium signaling, phagosome, cytokine-cytokine receptor interaction, folate biosynthesis, mucin type O-glycan biosynthesis, cell cycle, vascular smooth muscle contraction, progesterone-mediated oocyte maturation, and gap junction pathways were only significantly enriched by the DEGs in E26_vs._P2.
As shown in Figure 2D, the neuroactive ligand-receptor interaction pathway was commonly significantly enriched by the DEGs in each pairwise comparisons. However, the expression profiles of these DEGs enriched in this pathway were different, as shown in Figure 2E. Besides, some DEGs related to hormone synthesis and secretion was also significantly enriched in the neuroactive ligand-receptor interaction pathway, and their expression in geese pituitaries at four developmental stages were also verified by qRT-PCR. As shown in Figure 2F, expression of almost all these selected DEGs displayed changes in the same direction with those observed using RNA-seq, indicating the true reliability of our Illumina sequencing methods. Expression of melaninconcentrating hormone receptor-2 (MCHR2) gradually increased during the peri-hatch period. Levels of prolactin (PRL) mRNA increased from E22 to P2 but underwent a decline on P6. Expression of follicle stimulating hormone b (FSHb), thyroid stimulating hormone (TSHb), and thyroid stimulating hormone receptor (TSHR) genes decreased from E22 to P2 until an increased on P6.
Overview of the miRNA Transcriptome of Goose Pituitary During the Peri-hatch Period As shown in Supplementary Table S2, 66.88%7 4.73% of the miRNA sequencing reads were aligned to the goose reference genome. A total of 1,882 miRNAs were obtained from all samples, including 1,268 known-miRNAs and 614 novel-miRNAs. The length of known-miRNAs and novel-miRNAs were mainly distributed between 20 and 24nt, of which the miRNAs with a length of 22nt length were the most ( Figure 3A). The first base at the 5'end of the known and novel miRNAs were shown to have a strong bias toward 'U' ( Figure 3B). There were 78, 244, and 21 DE miRNAs identified between E22_vs._E26, E26_vs._P2, and P2_vs._P6, respectively. The E22_vs._E26 and P2_vs._P6 groups shared more than 50% of the DE miRNAs with the E26_vs._P2 group ( Figure 3C). By predicting the target genes for these DE miRNAs, 1,302, 3,030, and 2,037 target genes were found in groups E22_vs._E26, E26_vs._P2, and P2_vs._P6, respectively, and 418 target genes were overlapped in three groups ( Figure 3D).

Integrated Analysis of Differentially Expressed mRNAs and miRNAs Between Goose Pituitaries at Different Developmental Stages
These 418 shared predicted target genes of DE miR-NAs in three pairwise comparisons were significantly enriched in the ECM-receptor interaction, focal adhesion, and phosphatidylinositol signaling pathways (Supplementary Table S3). As shown in Figures 4A-4C, both the ECM-receptor interaction and focal adhesion pathways were commonly significantly enriched by the predicted target genes of DE miRNAs in each pairwise comparison. Furthermore, the interaction networks between 19 DE miRNAs and their 75 predicted target genes involved in the ECM-receptor interaction and focal adhesion pathways were constructed ( Figure 4D). Of them, three miRNAs including novel_miRNA_467, novel_miRNA_154, and novel_miRNA_340 had the most number of target genes, which might be hub miR-NAs. Moreover, the correlation between expression of DE miRNAs and differentially expressed target genes were also analyzed. Among them, expression of three miRNAs including novel_miR_120, tgu-miR-92-3p, and novel_miR_398 showed significant negative correlations with that of LAMC3 ( Figure 4E, Supplementary  Figure 2). In view of the central role of the ECM-receptor interaction pathway, we further verified the expression levels of the collagen and laminin gene family (COL2A1, COL4A4, COL4A5, COL4A6, COL6A6, and LAMC3) as well as 3 other DEGs (THBS2, THBS4, and VWF) involved in the ECM-receptor interaction pathway. As shown in Figure 4F, expression of almost all these selected DEGs determined by qRT-PCR showed changes in the same direction as those observed using RNA-Seq, further demonstrating the true reliability of the Illumina sequencing method. In addition, it should be noted that the expression profiles of the collagen gene family members were different. Among them, expression of COL6A6 increased from E26 to P2 but the others gradually decreased during the peri-hatch period.

DISCUSSION
Results from both PCA and heatmap analysis suggested that the pituitary transcriptome expression profiles were much more similar between E22 and E26 (during the embryonic period) as well as between P2 and P6 (during the post-hatch period). Meanwhile, it  PITUITARY TRANSCRIPTOME FROM EMBRYO TO HATCH was also observed that the number of both DEGs and DE miRNAs were the largest in E26_vs._P2 among 3 pairwise comparisons. These results altogether suggested that the goose pituitary transcriptome changes remarkably throughout the peri-hatch period, especially during the critical embryo-to-hatchling transition. Of note, 2 DEGs commonly present among 3 pairwise comparisons (Anser_cygnoides_newGene_32456 and LOC106031011) were differentially expressed throughout the peri-hatch period. As a newly identified gene, there currently lacks understanding of the functions of Anser_cygnoides_newGene_32456.
LOC106031011 that encodes the cystine/glutamate transporter-like protein can catalyze the transport of L-glutamic acid coupled with Na + and K + . Different types of glutamate receptors were shown to be present on the gland cells of the anterior pituitary (Petrusz, 1994;Villalobos et al., 1996). The latest studies found that glutamate mediated the apoptosis of prolactin-secreting cells and growth hormone-producing cells by activating type II metabolic glutamate receptors (Pampillo et al., 2002;Caruso et al., 2004). Therefore, LOC106031011 could be considered as one of the candidate genes regulating early goose pituitary development.
Based on the functional enrichment analysis of the DEGs in pairwise comparisons, it was proposed that the goose pituitary may help the individual adapt to sudden environmental changes during the peri-hatch period possibly by shaping its morphology and physiological functions, because the mostly-enriched pathways were related to tissue remodeling and regulation of cellular functions. Specifically, both the AGE-RAGE signaling and apelin signaling pathways were involved in vascular remodeling (Liu et al., 2016;Luo et al., 2018;Mughal and O'Rourke, 2018), but genes in these 2 pathways differ only in E22_vs._E26 and E26_vs. _P2 group. It was shown that AGE and its interactions with RAGE induced vascular hypertrophy by promoting the accumulation of ECM and stimulating the proliferation while inhibiting the apoptosis of vascular smooth muscle cells (Prasad, 2019). The versatile functions of the pituitary also depend on the coordinated development of the vascular networks (Scully et al., 2016). In addition, platelets play a key role in regulating of angiogenesis (Walsh et al., 2015). As the genes related to platelet reaction, including VWF, THBS2, and THBS4, their expression levels vary drastically from E22 to P2, while be stable between P2 and P6, which indicated that the development of goose pituitary blood vessels was before P2. Importantly, the function of goose pituitary seems to be drastically changed from E26 to P2. The embryonic to hatchling transition is marked by changes in the nutrient supply from yolk during the embryonic stage to feed after hatching, which is undoubtedly accompanied by the metabolic function changes, such as altered folate biosynthesis and mucin type O-glycan biosynthesis (Jing et al., 2013). Of note, the specific enrichment of the DEGs identified between E26 and P2 in the progesteronemediated oocyte maturation pathway indicated that the maturation of goose oocytes could be closely related to developmental changes in the pituitary during this period. Previous studies have shown that the thickness of ovarian cortex and the volume of the oocytes increase after hatching (Hu et al., 2020). These results altogether indicated that the remodeling of the goose pituitary may take place before P2, and its functional changes were particularly dramatic during the embryoto-hatchling transition.
In contrast to the KEGG pathways specifically enriched by the DEGs between different developmental stages, the neuroactive ligand-receptor interaction pathway was commonly significantly enriched by the DEGs in each pairwise comparison. Neuroactive substances are known to play a broad role in almost all physiology and organ development (N€ assel and Zandawala, 2019). These DEGs in pairwise comparisons were involved in different categories of neuroactive substances, including rhodopsin like amine, peptide, hormone protein, nucleotide like, secretin like, and metabotropic glutamate/pheromone. As the largest class of neuroactive substances, peptide hormone can affect the growth and development or reproduction of animals in a wide range of physiological activities (van den Pol, 2012;Nusbaum et al., 2017). In particular, most of genes related to hormone synthesis and secretion involved in the neuroactive ligand-receptor interaction pathway were differentially expressed in the goose pituitary during different developmental stages, including PRL, FSHb, TSHb, TSHR, and MCHR2. In poultry, PRL is generally accepted as the crucial regulator of the onset and maintenance of broodiness (Wilkanowska et al., 2014). Previous studies showed that a progressive expression of PRL in the pituitary occurs 1-2 days before hatching (Ishida et al., 1991). Meanwhile, the increase in circulating PRL concentration has been suggested to cause ovarian regression and broodiness (Porter et al., 1991;Jiang et al., 2005;Edens, 2011). These results suggested that changes in expression of PRL may be associated with physiological changes during the embryo-to-hatchling transition. In addition to PRL, expression of FSHb, TSHb, TSHR, and MCHR2 also varied in the goose pituitary during the peri-hatch period. Since the development of goose ovarian primary and secondary follicles begins after hatching (Hu et al., 2020), increased expression of FSHb in the pituitary from P2 to P6 indicated that FSH synthesis and secretion by the pituitary could be essential for early goose ovarian follicle development. Previous studies have shown that thyroid hormone is required for growth and development of poultry (King and May, 1984). Thyroid hormone may be the most important hormones in establishing a state of metabolic readiness to promote the occurrence of thermoregulatory responses in the early postnatal period (McNabb et al., 1984). However, the expression levels of TSHb and its receptor TSHR in different pituitary development stage were significantly different, and the expression trends of TSHb and TSHR during the peri-hatch period were similar. These data indicated that thyroid hormone may act on the pituitary gland through negative feedback after acting on other target tissues during the embryonic to posthatching period. Although MCH has not been annotated in the genome of goose (Lu et al., 2015), their receptor (MCHR2) has been detected to be differentially expressed during the peri-hatch period. MCHR2 is considered to be a functionally redundant pseudogene in chickens (Cui et al., 2017), its differential expression during the peri-hatch period seems to indicate its function in goose pituitary development. These results indicated that the function of the pituitary has undergone significant changes during early developmental stages, which are strongly associated with differential expression of genes involved in the neuroactive ligand-receptor interaction as well as hormone synthesis and secretion pathways.
Both the target genes of DE miRNAs and the 418 common target genes in pairwise comparisons were significantly enriched in the ECM-receptor interaction pathway. Previous reports have shown that ECM proteins are important for the normal morphology development in the pituitary (Paez-Pereda et al., 2005). ECM components can affect the function of FS cells and play an important role in guiding local cell arrangement in the anterior pituitary (Horiguchi et al., 2010). The individual components of the ECM, such as laminin, fibronectin and collagen, have effects on cell proliferation, differentiation, morphogenesis, and hormone production (Lukashev and Werb, 1998). Expression of multiple collagen gene family members, including COL2A1, COL4A4, COL4A5, COL4A6, and COL6A6, was significantly during different during the peri-hatch period. As a component of ECM, laminin maintains tissue structure and transmits biological information for cells (Aumailley and Gayraud, 1998). In this study, expression of LAMC3 and LAMA4 (laminin subunit alpha 4) were significantly different during the embryo-to-hatchling transition, but their expression profiles were opposite. Previous studies have shown that LAMC3 is important for the development of the human brain (Barak et al., 2011), but there is no research on the role of LAMA4 in the brain or pituitary development. Therefore, whether LAMA4 can compensate for the function of LAMC3 in the goose pituitary after hatching remains to be further studied. It is known that miRNAs usually inhibit the expression of their target genes. LAMC3 has a significant negative correlation with the expression of novel_miR_120, tgu-miR-92-3p, and novel_miR_398, which has a predicted targeting relationship. At the same time, the expression of novel_miR_120, tgu-miR-92-3p, and novel_miR_398 were relatively high in miR-NAs, these results indicated that novel_miR_120, tgu-miR-92-3p, and novel_miR_398 may be play a key role in goose pituitary development by targeting LAMC3. These results suggested that these identified DE miR-NAs regulate the pituitary remodeling and functions possibly by targeting genes involved in the ECM-receptor interaction pathway.
In conclusion, this study represented the first to systematically describe the mRNA and miRNA transcriptome profiles in goose pituitary development during the peri-hatch period. These identified DE mRNAs and miRNAs as well as their predicted interaction networks involved in regulation of tissue remodeling and cellular functions were most likely to play critical roles in facilitating the embryo-to-hatchling transition. These results provide novel insights into the early developmental process of avian pituitary and will help better understand the underlying molecular mechanisms.