Network Biology led identication of critical miRNA modules for High Altitude hypoxia and their potential as dietary supplements

Early ascent to high altitude can cause severe damage to body functions and may lead to many fatal high-altitude disorders. To cope up with such conditions, the human body undergoes physiological and biochemical changes in order to adapt to extreme environmental conditions at high altitudes. Several microRNAs (miRNAs), Transcription Factors (TFs), and genes have been studied separately for their role in early adaptive molecular responses. We hypothesize that network analysis of miRNA-TF-gene coregulatory networks of circulatory miRNAs (CmiRNAs) which are differentially expressed at high altitudes could reveal a complex regulatory functional module that might be controlling molecular adaptive responses at high altitude. A comprehensive and non-redundant list of differentially expressed human CmiRNAs during high altitude ascent was collated and 470 Feed-Forward Loops (FFLs) tripartite motifs were identied in the miRNA-TF-gene co-regulatory networks. Network analysis and K-means clustering identied 11 biologically overrepresented FFLs regulated by 8 miRNAs hsa-miR-335-5p, hsa-miR-26a-1-3p, hsa-miR-210-3p, hsa-miR-193b-3p, hsa-miR-17-5p, hsa-miR-16-5p, hsa-miR-5582-5p and hsa-miR-130a-3p. Pathway enrichment identied metabolism and inammation as important hallmark responses responsible for high altitude adaptation. These miRNAs were evaluated for their supplementation from dietary sources. Phylogenetic analysis with 4 other mammalian and non-mammalian species showed that Bos taurus (cattle) milk could be a possible source of these dietary miRNAs. Exogenous miRNAs bta-mir-16-1, bta-mir-130a, bta-miR-335, and bta-miR-210 have >95% of sequence similarity and could become potential dietary miRNAs candidates. The sequence, structural properties, and high AGO2 binding eciency of all these exogenous miRNAs show good serum stability and cellular uptake possibility in the mammalian host. Bta-miR-210 with highest Minimal Folding free Energy Index (MFEI) and highest Atomic Contact Energy (TF)-PPARG (target gene). Eighth and ninth FFLs are hsa-miR-210-3p (miRNA) - HIF1A (TF)-BNIP3 (target gene) and hsa-miR-210-3p(miRNA)-HIF1A(TF)-LDHA (target gene) respectively as they contain one common exclusive interaction hsa-miR-210-3p-HIF1A. Hence, all interactions associated with these miRNAs were subjected from high altitude specic miRNA-TF gene coregulatory network. The TF and genes targets of this module were functionally annotated to identify key molecular mechanism regulated by the four miRNAs.

these physiological and biochemical changes, emphasizing the contribution of the underlying complex molecular system in high-altitude adaptive phenotypes [11,12]. TFs play an important role in regulating biomolecular processes and gene expression during high altitude hypoxia adaptation [13]. HIF (Hypoxia Inducible Factor) is the most important TF that plays a major role in the activation of a multitude of hypoxia stress responses [13]. Recent research further indicates that polymorphism in hypoxia, nuclear receptor transcription factors i.e. Endothelial PAS domain-containing protein 1 (EPAS1), Egl-9 Family Hypoxia Inducible Factor 1 (EGLN1), and Peroxisome proliferator-activated receptor gamma (PPRAG) have also been proposed to be associated with acclimatization of an individual [11]. Recent reports have also found altered miRNAs expression in humans during acute ascent to high altitude and this change in miRNAs expression has been correlated with adaptive molecular response [14][15][16]. miRNAs are small noncoding RNAs that are known to regulate gene expression post-transcriptionally and translationally in various physiological and pathophysiological conditions [17]. They have been found to be stable in different body uids like blood serum, plasma, breast milk, lymph uid [18]. Some of the miRNAs are secreted into the extracellular space and then transported to the circulating body uid like peripheral blood known as cmiRNAs [19]. Although, these reports have identi ed altered miRNA pro les in both acclimatized and acclimated individuals at similar altitude and environmental conditions. But all these altered cmiRNAs pro les are different from each other and govern different molecular functions. The altered miRNAs expression, physiological and biochemical traits contributing to early acclimatization are strongly correlated [14,20]. Thereby opening an opportunity to identify the underlying complex regulatory mechanisms, governing molecular responses for high altitude adaptations.
Hence, we designed a meta-study to collectively review these pro les and elucidate their common regulating molecular mechanism during high altitude ascent. The study has intended to identify a set of differentially expressed CmiRNAs that are having a high correlation with the high-altitude adaptation molecular mechanism. The comprehensive list of 299 CmiRNAs was curated from these literature studies and their experimentally-validated miRNA-target interactions (MTI) were identi ed from validated repositories. Tripartite miRNA-TF-gene co-regulatory networks were built and analysed by identifying the 470 FFLs. FFLs are interconnecting patterns that repeat in many different parts of a network and form key functional modules [21]. FFLs have been demonstrated as one of the most important motif patterns in gene regulatory networks that govern various aspects of normal cell functions and diseases [22].
Detailed examination of the high-altitude studies shows that 6 of the nine key miRNAs are relatively under expressed in low landers as compared to the high altitude native population. Evidence shows that cmiRNAs can be compensated with exogenous transfer miRNAs that are derived from dietary substances [24][25][26]. These dietary miRNAs could be absorbed through the intestine and remain active to regulate gene expression using the host machinery [25,26]. To identify potential dietary miRNA candidates for compensating 6 under-expressed cmiRNAs, phylogenetic analysis with other mammalian and nonmammalian species was performed. The identi ed exogenous source miRNAs have checked for their cross-kingdom serum stability and cellular uptake for the potential to be supplemented as dietary miRNAs to compensate for the downregulated complementary human miRNAs. For this, sequence and structural properties of these exogenous miRNAs have been compared with gold standard dietary miRNA OSA-MIR168a for cross-kingdom serum stability and miRNA-AGO2 docking for cellular uptake. Figure 1 provides a graphical view of the methodology which is followed in construction and analysis the miRNA-TF-gene coregulatory network to identify biologically overrepresented functional modules. Methodology 1. Collection of differentially expressed miRNA at high altitude.
To accumulate a comprehensive list of miRNA, data mining was performed by using PubMed clinical queries with keywords "high altitude" and "miRNAs". This query was run in ve different categories i.e. 'therapy', 'etiology', 'diagnosis', 'prognosis' and 'clinical prediction guides' with 'broad' scope to imbibe a larger amount of studies [27]. The query yielded an exhaustive list of PMIDs having differentially expressed miRNAs at high altitude. The list was further clustered together and duplicity was removed to get a non-redundant test list of PMIDs for ranking the text. Data mining tool, Medline Advance Ranker was used to rank the list of PMIDs based on the training set parameters [28]. The weighted discriminative words become the scoring function of the training model. The model was run against the test set, which further ranked the PMIDs according to the scoring function described in the training model. The data was plotted in R (PubMed Ids vs |log (P)|) with abline.
2. Identi cation of miRNA-gene target and miRNA-TF target pairs.
The miRNA-target gene pairs were identi ed using the MultimiR R package [29]. The package collects the data from several validated databases, including mirTarBase [30], TarBase, and miRecords [31].To increase the stringency, only highly validated targets (i.e. validation from reporter assay, luciferase assay, proteomics, western blot, and qPCR) were retained in the study. All human TFs and TcoFs set was formed by combining data from different databases i.e. TFcheckpoint [32], DBD [33], ORFeome [34], TcoF-DB V2 [35], TFCat [36], TFClass [37], TRANSFAC [38]. Further, this all human TF and TcoF set was used to label miRNA targets. Those matching the datasets were labelled as TF and those that did not match were labeled as the genes.
3. Construction of miRNA-TF co regulatory network.
After deciphering the miRNA targets into gene and transcription factors to identify miRNA-gene and miRNA-TF interactions, other interactions such as TF-gene and TF-miRNA were also added from respective public repositories. The TF: TG interactions were fetched from OregAnno 3.0 [39] and TRRUST V2 [40] databases. TF: miRNA interactions were fetched from TransmiR [41] and PuTmiR [42] databases.
The nal le was compiled into two columns (.SIF) format and was used to visualize the network in Cytoscape [43].

Functional Annotation.
Database for Annotation, Visualization and Integrated Discovery (DAVID) is a powerful tool to mine functions of gene target [44]. The target miRNA-TF-gene regulatory network was submitted to DAVID to identify their Gene Ontologies (GOs). GO functional enrichment analysis yielded a large number of GO terms. As there are a large number of GO sets and most of them are completely overlapping to other smaller GO sets. Hence to analyse them together, Gene set enrichment analysis (GSEA) was performed using the 'Enrichment map' plugin of Cytoscape [45]. The GSEA yielded enriched GO terms and the result was shown using heatmap. Further, the enriched GO terms were distributed using boxplot to identify exclusive GO sets that range from exclusive to highly exclusive. 5. 5. Identi cation of miRNA -TF co-regulatory FFLs motifs.
The miRNA-TF-gene coregulatory network was analysed by identifying FFL motifs from the network. The FFL network module was identi ed using the 'motif discovery' plugin in cytoscape and the tight module Clique Percolation Method (CPM) of C nder tool [46]. Randomization test was performed to evaluate the signi cance of the FFLs observed in the set of TFs, miRNAs and genes. For this, we compared how often one FFL appears in the real network to the number of times it appears in randomly generated networks formed by degree preserving randomization algorithm. In order to retain biological key driver nodes, we used a degree preserving randomization algorithm of the 'igraph' R-package. Randomization process repeated 1000 times. Z score for each motif type was calculated to identify the signi cant motifs.

Zscore=No−Nm/σ
Here, No is the number of motifs observed in the real network, whereas Nm, and σ are the mean and standard deviation of the motif occurrence in 1000 random networks, respectively. Five different FFL subgraph (motifs) were identi ed in the miRNA-TF-gene coregulatory network 7. Motifs evaluation.
The FFL network module was evaluated on the basis of two Biological Network Motifs (BNM) algorithm i.e. EDGEGO-BNM and EDGEBETWEENNESS-BNM algorithms [23].
EDGEBETWEENESS-BNM-This algorithm uses edgebetweeness topological parameters to identify biologically signi cant edges of the FFL motifs. The edgebetweeness centrality of each edge in the network was calculated. Edgebetweeness centrality de nes the number of shortest paths that passes through a vertex. This topological parameter helps in identi cation of the bottleneck interaction.
EDGEGO-BNM-In this algorithm, miRNA target genes and TFs of the FFL network were rst functionally enriched with gene ontology to identify biological relevant network motifs. We have modi ed the algorithm by performing pathway enrichment analysis instead of gene ontology because pathways enrichment analysis yields more signi cant terms than gene ontology. For this, the two target sets of miRNA-TF co-regulatory networks i.e. genes and TFs were separately used for enrichment. The Reactome pathway database was used for pathway enrichment. Quartile distribution was used to identify any possible outliers in both pathway sets. For this, the mode of log [p.values] obtained after pathway enrichment was used to draw two boxplot of both pathways sets. The gene set of the most outlier and its interaction miRNA, TF, gene partners were used to construct EDGEGO-BNM network. Lastly, the edge between each edge in the EDGEGO-BNM network was calculated to identify bottleneck interaction from the network.
The edgebetweeness centrality is an edge related topological parameter. There are thirteen other node related topological parameters by which the EDGEGO-BNM network can be assessed. For this, K-means clustering was performed for all nodes of EDGEGO BNM network using their normalized topological parameters [47]. K was determined to know the number clusters in K-means clustering by using three different methods i.e. Elbow method: 4 clusters solution suggested, Silhouette method suggested 6 clusters and Gap statistic method only one cluster. As all three-methods predicted three different results, So NbClust R package was used to compute frequency of clusters in 30 different methods simultaneously [48]. Before performing the K-means clustering, variable contribution in dimension 1 and 2 was determined to identify major variables contributing in the clustering. Factoextra and FactoMineR R Packages were used for performing K-means clustering and ggplot2 R package was used for its visualization [49]. Further, gene ontology enrichment of TF/gene present in these clusters are being performed.
8. Prevalence of molecular responses during acclimatization at high altitude as hallmark.
Prevalence of molecular responses during acclimatization at high altitude was calculated using interactive learning approach to classify them as hallmark. Data mining was used for tagging literature with appropriate molecular responses. To develop a training model for classi cation of the molecular responses as hallmark, a training dataset was constructed. The training dataset comprised of 41 literature PMIDs that were mined using keyword "acclimatization at high altitude", "adaptation at high altitude". Unsupervised learning approach was used to train the model and the scoring function was developed based on the discriminative word's weightage [28]. The model was validated using Leave-oneout cross-validation algorithm and predicted performance parameters were calculated during validation [50]. The receiver operating characteristics curve (ROC) shows that model is highly sensitive. The test datasets were developed again by using molecular pathways regulated by 9 key miRNAs in high altitude speci c miRNA-TF-gene coregulatory network as keywords. This will help to check if these molecular responses are apt candidate for the hallmark responses with negative and positive control. The 14-test dataset consisted of PMIDs mined using keywords 'body weight homeostasis', 'metabolism of angiotensinogen to angiotensin', 'erythropoiesis', 'cell development and senescence', 'circadian rhythm', 'hypobaric hypoxia regulation', 'nuclear receptor factor', 'HIF-1 responsive genes expression', 'cytokines signaling', 'cell cycle regulation', 'adherence and junction genes', 'protection during oxidative stress', 'positive control' and malarial as negative control'. These datasets were tested against trained model to classify them as hallmark response during acclimatization at high altitude and scored Log (AUC) was plotted in R.

Phylogenetic analysis
To perform phylogenetic analysis, a comprehensive list of differentially expressed human (Homo sapiens), Tibetan pig (Sus scrofa), Cattle (Bos tauras) Tibetan yak (Bos mutus), chicken (Gallus gallus) miRNAs at high altitude was obtained from published studies [51,52]. Individual lists of miRNAs from the above mentioned species were made to intersect to obtain a nal list of 100 common miRNAs. To perform phylogenetic analysis, precursor sequences of 100 miRNAs were retrieved, aligned, analysed and clustered to establish a phylogenetic relationship between them. 9.1 Sequence retrieval-The source for all microRNA sequences in the One hundred precursor sequences of human miRNA, were retrieved from miRBase [53], mirGenedb [54] and dbDEMC 2.0 [55] database in Fasta format.
9.2 Sequence alignment-BLAST (Basic Local Alignment Search Tool) [56] was performed for each of the human (Homo sapiens) miRNA sequences retrieved from NCBI using the online Geneious 11.0.5 software [57]. Precursor sequences were studied for their similarity patterns by feeding the FASTA format into the software. Each of the human miRNA sequences was studied for their similarity patterns and BLAST was therefore performed by feeding the FASTA format of sequence.
9.3 Phylogenetic analysis-Phylogenetic analysis of each human miRNA sequence through UPGMA (Unweighted Pair Group Method with Arithmetic Mean) was performed using Geneious 11.0.5 software [57]. A tree was constructed by the software showing the ancestral relationship among the sequences. The tree gives different clusters, showing their relationship with each other. The sequences which lie in the same cluster are more closely related than those lying in distant clusters. 9.4 Hierarchical Clustering-To expand the results of phylogenetic analysis for better understanding hierarchical clustering was performed on data set. To perform hierarchical clustering the 'heatplus' R package was used. The results were visualized in the form of heatmap.

Statistical Analysis
All statistical analyses were performed with the R software (version 3.2.3, R Foundation for Statistical Computing, Vienna, Austria). Package 'ggfortify' / 'ggplot2' used visualization of plots [58]. 'Boxplot' R package was used for Quartile distribution [59]. 'Heatplus' R package for constructing heat map 11. Molecular docking between miRNA-and AGO protein Argonaute (AGO) protein is a family of protein which assists miRNAs to bind with mRNAs of the target genes. Hence, study of the binding mechanism between AGO protein and miRNAs is important for understanding their RNA silencing mechanism. For this, secondary structures of proposed exogenous miRNAs were predicted using RNAfold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) to identify dot bracketed structures of miRNA [60]. The resulted dot bracket structures were further used for prediction of their tertiary structure using RNA COMPOSER (http://rnacomposer.cs.put poznan.pl/) [61]. Further, the three-dimensional structure of AGO protein was retrieved from Protein Data Bank (PDB ID: 3F73) and was prepared for docking through energy minimization and water molecules removal though Discovery Studio 3.5 suite (http://accelrys.com/products/discovery-studio/visualizatation-download.php) [62]. Docking between three dimensional structures of protein and ligand, may be a biomolecule like miRNA or a gene is an e cient computational method to inspect the molecular interaction. Docking were performed between AGO protein and miRNA (bta-mir-16-1, bta-mir-130a, bta-miR-335, bta-mir-210, positive and negative control) using PatchDock web server (http://bioinfo3d.cs.tau.ac.il/PatchDock/), algorithm ranked the docked complexes on the basis of highest geometrical shape complementary binding energy scores [63].

miRNA cross-kingdom stability
Recent studies have formulated criteria for identi cation of potential cross-kingdom transferred miRNAs that are stable in human serum by taking miR-168 as gold standard such as a) miRNA secondary structures were selected on the basis of higher minimal folding free energy (MFE) and MFEI, MFEI usually being over 0.85 [64,65]. b) Structures with not more than 9 mismatches, maximum 4 bulges and no more than 3 loops were considered [65]. Mfold server was used to identify all these properties of the four sequences i.e. bta-mir-130a, bta-mir-16b, bta-mir-210, bta-mir-335.

Results And Discussion
A network-based approach is used to collectively study differentially expressed human cmiRNAs during high altitude ascent. As mentioned in the methodology, the literature associated with high altitude cmiRNAs was extracted from the PubMed database. This literature dataset was ranked based on the training model that was developed using discriminative word clouds of high-altitude adaptation. The 299 cmiRNAs are curated from the four [14,16,20,67] out of six PMIDs and. 959 human volunteer samples were analysed, which includes native Tibetans -362, newly arrived migrants to high altitude -126, native low landers -233 and 238 long-time migrants to high altitude.
These 299 differentially expressed high altitude speci c human CmiRNAs experimentally validated target interactions were identi ed using "MultimiR 2.0" R package [29]. The miRNA-target interactions were further classi ed as miRNA-gene and miRNA-TF interaction by labelling miRNA targets as TFs or genes using all TFs and TcoFs dataset. The nal datasets contain 3455 miRNA targets with 2991 genes and 464 TFs respectively. In terms of interactions, there are 8519 miRNA-genes and 1655 miRNA -TF interactions. To construct miRNA-TF-gene co-regulatory network 364 TF-miRNA and 566 TF-gene interactions were added from their validated databases as mentioned in methodology. The nal directed tripartite miRNA -TF gene co-regulatory network contains 3754 nodes (299 CmiRNAs, 509 TFs and 2936 genes) and 13407 interactions and visualized using Cytoscape 3.2.1 [43]. Degree distribution topological analysis was performed on the miRNA-TF gene coregulatory network showed that it follows the power law and network has scale free topology. The scale-free topology is very common in biological networks.

Functional Annotation and Gene Set Enrichment Analysis
The targets of miRNA-TF-gene coregulatory network was functionally annotated with Gene Ontology (GO) using DAVID server. The functional enrichment analysis of 2991 target genes in the network produced 3300 GO terms. These terms were again subjected to Gene Set Enrichment Analysis (GSEA), yielded 463 enriched sets/clusters. Further quartile distribution based on enrichment score of 463 enriched sets has yielded many outliers ( Figure 2). The outliers can be divided into two categories based on their grouping exclusive terms and highly exclusive terms. The exclusive terms can be seen as the extended part of data in the fourth quartile ( Figure 2). But highly exclusive terms were so far from mainstream quartiles that it can't be included. Hence, these highly exclusive terms can be represented as most enriched gene sets in the network. The 5 most enriched sets that were obtained are phosphoprotein/protein kinase, cellcell adhesion junction proteins, membrane proteins, glycoproteins, tyrosine protein kinases. These topmost terms are phosphoproteins/protein kinase with an enrichment score of 61.33. This GO term has been apparently observed in many GO analysis of other hypoxia datasets because protein kinases are major signal transducing molecules that can transduce hypoxic stress signals in the cytoplasm to the nucleus, where activation of many hypoxia responsive factors such as HIF1, NFKB takes place. Other clusters that appeared were cell adherence proteins and membrane proteins with second and third best enrichment score of 38.7 and 28.7 respectively. FFL analysis.
Biological systems are often built of recurring circuit modules that carry out key functions. Biomolecules such as miRNA, transcription factors and genes follow these recurring patterns in co-regulatory networks that regulate the molecular and cellular responses of living cells. These recurring patterns are called motifs, FFL is the most common type of motifs that are found overrepresented in most biological networks especially in regulatory networks. These FFL motifs consist of three entities miRNA, TF and gene. The TF and its target gene are directly regulated by the miRNA, and simultaneously miRNA or TF regulate each other. Recently these FFLs have helped the researcher in identifying important regulatory mechanisms during many disease pathologies. Since miRNAs signalling involves the functioning of the enormously complex system of molecular interaction at the system level. Hence, identifying FFL motifs in complex networks could further reveal the regulatory modules that can be used to de ne particular ecological traits at high altitudes. FFLs were identi ed in the directed miRNA-TF gene co-regulatory network by using the 'motif discovery' plugin of Cytoscape. The network contains 470 FFLs, having 494 nodes (82 miRNA, 85 TF, 327 genes) and 2942 edges between them.

EDGEGO-BNM
FFL motif network is evaluated second time with EDGEGO-BNM evaluation. For this, rst pathways enrichment of genes and TFs in the FFL network have been performed that shows 1067 pathways and 27 pathways respectively. The quartile distribution of 1067 pathways terms from pathways enrichment of genes yielded with no outliers. Whereas, quartile distribution from pathways enrichment of TFs generated four outliers i.e. 'nuclear receptor transcription pathway' is the most signi cant term obtained in pathway enrichment of TFs followed by 'circadian gene expression', 'clock genes' and 'PPARA activates gene expression'. An outlier is a data point that differs signi cantly from other observations and in general, the enrichment p.values have been calculated by comparing the observed outcomes of an annotation term with the expected outcomes; individual terms that are beyond some limit or deviated from the original dataset said to be enriched. So, 'nuclear receptor transcription pathway' was adjudged as the most enriched pathway term and its gene sets (PPARA, RARG, PPARG, NR3C2, NR4A2, ESR1, RARA, AR and NR3C1) with their interacting miRNA, TF and gene partners were used to construct EDGEGO-BNM network. Nuclear receptors and their induced transcription factors are a consortium of proteins that are responsible for sensing hormones like steroid, thyroid and certain other molecules to regulate the speci c gene expression in the nucleus, thereby controlling the development, homeostasis, and metabolism of the organism. Nuclear receptors are found essential for maintaining body homeostasis and metabolism at high altitude [70,71] . The EDGEGO-BNM network has 23-miRNAs, 29-TFs, and 23-genes with 302 interactions (Figure 4a). Subsequently, edgebetweeness centrality of all its 302 interactions have been identi ed and on quartile distribution 7 outliers are obtained (Figure 4b). These are exclusive interactions in the network with hsa-miR-5582-5p -SP1, ESR1 -PTPRG, ESR1 -SP1, hsa-miR-335-5p -SP1, ERG -ESR1, ESR1 -TSC22D3, ESR1 -CCND1, hsa-miR-130a-3p -ESR1, hsa-miR-335-5p -NR4A2 are signi cantly enriched and bottleneck edges of the network.
Prevalence of molecular responses during acclimatization at high altitude as hallmark.
The molecular responses obtained in miRNA-TF-gene network module which is regulated by the differentially expressed cmiRNAs; hsa-miR-335-5p, hsa-miR-26a-1-3p, hsa-miR-210-3p, hsa-miR-193b-3p, hsa-miR-17-5p, hsa-miR-16-5p, hsa-miR-5582-5P, hsa-miR-130a-3p have been tested to be hallmark responses using knowledgebase ranking of biomedical literature. For this, the training dataset for classi cation of these molecular responses as a hallmark was built. The training dataset comprises 41 PMIDs that are mined using keyword "acclimatization at high altitude", "adaptation at high altitude" from PubMed database. The model was trained using Medline ranker advance to identify scoring function based on discriminative word weightage cloud. After performing leave one out cross-validation the model achieved 95.3% of accuracy. The model was used to rank the 13 molecular responses regulated by four cmiRNAs as hallmarks for high altitude adaptation. The ranking pattern of these molecular responses was compared with positive (training dataset) and negative (malaria dataset) control. Quartile distribution of the normalized p.values of ranked molecular responses related literature shows that 'cytokine signaling', 'interleukin signaling', 'cell cycle regulation', 'cell development and senescence', 'body weight homeostasis' and 'circadian rhythm' are normally distributed and their distribution is very close to the positive control. Also, recent reports have proved that these responses appear immediately after high altitude ascent. 'Erythropoiesis', 'nuclear receptor factor regulation', 'protection from oxidative stress' are slightly distant from positive control but slight variation could be possible due to larger sample size of the literature than the training set. Hypoxia responsive genes and adherence and junction genes have unequal data distribution in quartiles. Hence most of the data trend is towards quartiles that are closer to the positive control including outliers that are closest with the positive control trend ( Figure 5). Further comparative analysis Area Under Curve (AUC) showed a high correlation between 13 molecular responses and high-altitude adaptation. Hence, above results prove that these molecular responses are hallmark responses for high altitude adaptation.

Phylogenetic analysis
Phylogenetical analysis was performed between 102 common precursor-miRNAs sequences of human and some high-altitude living mammals including Sus scrofa, Bos mutus, Bos taurus and non-mammals Gallus gallus. Further hierarchical clustering of the percentage sequence similarity gave 16 clusters ( Figure 6). The sequences having the highest percentage of similarity are more closely related, while the ones with <90% similarity shows the dissimilarity among them. The cluster number 12 contains 6 human precursor miRNAs hsa-miR-497, hsa-miR-335, hsa-miR-212, hsa-miR-190b, hsa-miR-148b, hsa-miR-125b having high percentage similarity with all three mammals. Similarly, cluster number 11 contains 17 precursor miRNAs including hsa-miR-130a that have more than 93% of sequence similarity with mammals and <90% with Gallus gallus. Cluster 13 contains 6 miRNAs including hsa-miR-16 having >90% sequence similarity with all three mammals and Gallus gallus. Whereas, cluster 4 contains 5 miRNAs including hsa-miR-210 and shares >93% sequence similarity with Bos tauras, Bos mutus, Gallus gallus respectively. Lastly, cluster 1 contains 13 miRNAs that have sequence similarity <90% with both mammalian and Gallus gallus. Overall phylogenetical analysis have shown that human cmiRNAs hsa-miR-335-5p, hsa-miR-210-3p, hsa-miR-16-5p and hsa-miR-130a-3p identi ed in our network analysis have shared > 95% of sequence similarity with the bta-miR-335, bta-mir-210, bta-mir-16-1 and bta-mir-130a of Bos tauras species. These miRNAs are also found differentially expressed in bovine milk at high altitude. Hence, phylogenetic analysis supports our hypothesis that these cmiRNAs are part of regulating modules controlling early adaptive mechanisms in humans during high altitude ascent.
miRNAs cross-kingdom serum stability was checked by comparing them with sequence and secondary structure properties of gold standard miR-168a showed all proposed miRNAs have ≥0.85 MFEI score, <9 mismatches and have only 1 bulge in pre-sequence (Table 2). But only bta-mir-16b, bta-mir-210 and btamir-335 have ≤ 3 Interior loops, whereas bta-mir-130a precursor sequence contains 4 Interior loops ( Table   2). This makes bta-mir-16b, bta-mir-210 and bta-mir-335 as more potent exogenous miRNAs to have cross-kingdom stability. The exosomes and AGO2 binding are the most common mechanism studied for transfer and uptake of exogenous miRNAs from peripheral blood to cells. Docking studies of potential exogenous miRNAs are performed with human AGO2 protein showed bta-mir-210 with highest ACE (-999) score ( in bta-miR-193b and bta-miR-335 with aromatic rings were also found in the binding pocket of miRNAs within a distance of 3.0 Å.

Discussion And Conclusions
The prevalence study reveals 'cytokine signaling' and 'bodyweight homeostasis' as the two key molecular responses regulated by the module of 8 miRNAs with associated TFs and target genes have comparatively high occurrences in most of the high altitude studies than other responses. So, 'cytokine signalling' and 'bodyweight homeostasis' are closest to be the hallmark responses responsible for high altitude adaptation. The cytokines can work as both pro-in ammatory and anti-in ammatory factors to regulate the in ammatory response and in ammation-based immune response can be bene cial in reducing cellular and tissue damage due to oxidative stress and ischemic condition [72]. Whereas, in order to maintain 'body weight homeostasis' in high-energy demanding environments like high altitude leads to an increase in blood pressure for maintaining oxygen saturation in tissues to meet the energy demand [73]. Energy demand increases due to the initiation of a negative energy balance mechanism because of acute exposure at high-altitude, which is associated with the loss of appetite (Anorexia) and leads to weight loss [74]. Though, to understand the framework of bodyweight homeostasis various models have been proposed and symmetric or asymmetric biological control of body weight is the most common mechanism [75]. The most common of this model is adipokines leptin secreted from adipose tissue which helps in maintaining whole-body homeostasis [75]. Adipokines are the in ammatory factors produced and secreted by adipocytes e.g. Leptin, Adiponectin, Interleukin-6 (IL-6), Plasminogen Activator Inhibitor-1 (PAI-1), Tumor Necrosis Factor-Alpha (TNFα), Interleukin 8 (IL-8), Interleukin 10 (IL-10), Interleukin 1B (IL-1β) Interferon-gamma (IFN-γ), etc. Adipokines secretion takes place due to the activation/polarization of adipose tissue macrophages that causes an induced proin ammatory state [76]. A low grade in ammation called 'metain ammation' causing polarization of proin ammatory M1 macrophages and anti-in ammatory M2 macrophages have been linked to the development of adipose tissue in ammation. Macrophage polarization to M1 state can lead to chronic in ammation causing metabolic syndrome like obesity [76]. Literature reports show anorexia, weight loss, and alterations in metabolism at high altitudes have been associated with adipokines release e.g. Leptin release is closely correlated to anorexia and weight loss at high altitudes [77,78]. Similarly, adipokines like acylationstimulating protein (ASP), adiponectin and IL-6 have been proposed as a biomarker for alterations in metabolism due to high-altitude ascent at 5300m [78]. Also, the relationship between excess adiposity and increased blood pressure is well established [73,79]. The meta-in ammation process has been widely studied in obese individuals and obesity is a well-established risk factor for high-altitude diseases like AMS), especially for the Tibetan plateau region [80,81]. This shows that in ammatory and metabolic mechanisms are somehow working together for high altitude responses and the FFL module of 8 miRNAs and associated TF and genes are part of the regulatory nexus controlling this mechanism through adipokines.
The cross-talk between miRNAs and adipokines has been known to mediate the in ammatory and metabolic state of the body to maintain core internal temperature, blood pressure, and body weight homeostasis [82]. A recently published report on adipose secreted miRNA miR-210 shows that it can govern important high-altitude hypoxia pathways. It inhibits adipose browning to reduce the over secretion of pro-in ammatory cytokines and help to curb the occurrence of diseases like obesity, diabetes, and other cardiovascular disorders [83]. Hsa-miR-210 has appeared as one of the biologically overrepresented miRNAs during our analysis. Similar association proofs of other 8 miRNAs with adipokines secretion and in ammatory responses that have been identi ed during FFL analysis are also available in the literature. Some of these are; hsa-miR-16 and hsa-miR-335 also known as adipocytespeci c miRNA. Hsa-miR-16 is associated with a change in body weight whereas, miR-335 is involved in adipose tissue in ammation [84,85]. MiR-17 controls adipocyte expansion and differentiation that resulted in a shift toward the dominance of pro-in ammatory adipokines released over anti-in ammatory adipokines [86]. Similarly, miR-193b and miR-130 can also regulate the expression and the release of adipokines [87,88]. Hsa-miR-26a-1 has been identi ed as an inhibitor of a part of adipocyte differentiation and a novel therapeutic candidate for obesity [89]. Lastly, no association of hsa-miR-5582-3P with adipokines have been reported to date. But it has been found to regulate other proin ammatory cytokine receptors like TGFβ-R1 and TGFβ-R2, to modulate TGFβ signaling [90].
Regulation of Adipokines release to manage the in ammatory and metabolic state of the body is a complex mechanism and found to be under the strict control of miRNAs and TFs during high altitude stress. The biologically overrepresented FFLs identi ed in FFL network analysis can be used to represent these complex regulatory systems because all these FFLs contain at least 2 of its three interactions as bottleneck interaction. In 4 of the 11 biologically overrepresented FFLs i.e. hsa-miR-16-5p -SP1-VEGFA, hsa-miR-335-5p-SP1-APP, hsa-miR-5582-5p-SP1-NR3C1, hsa-miR-335-5p-EGR1-SP1, SP1 has been a coregulating partner of hsa-miR-5582-5p and hsa-miR-335-5p regulation 3 target genes Vascular endothelial growth factor (VEGFA), Amyloid Precursor Protein (APP), Nuclear Receptor Subfamily 3 Group C Member 1(NR3C1), and regulated by TF Early Growth Response 1 (EGR1). Moreover, the FFL network evaluation through EDGEGO-BNM shows that SP1 is there in 3 out of the ve bottleneck interactions.
SP1, a ubiquitously-expressed transcription factor that plays a vital role in the regulation of numerous genes required for normal cell functions like cellular processes, including cell differentiation, cell growth, apoptosis, immune responses, response to DNA damage, etc [91]. SP1 can directly or indirectly regulate the release of different adipokines e.g. IL-1, CCL2, RBP4, VEGFA, SERPINE4, TNF-α, etc. and control cytokine implicated in in ammation, apoptosis, and cell survival as well as induction of insulin resistance [92]. SP1-APP is one of the four essential bottleneck interactions in the FFL network. APP is a microglial acute-phase protein. These acute-phase proteins are established markers for an in ammatory response [93]. Many studies have projected a high correlation of acute-phase proteins increase with in ammatory cytokines or adipokines such as IL-1β, IL-6, and TNF-α. APPs are potentially known to be biomarkers for the diagnosis of AMS [93]. SP1-VEGFA is the second bottleneck interaction, VEGFA is mainly stimulated by hypoxia and is an important signaling protein involved in hypoxia-stimulated vasculogenesis and angiogenesis. VEGFA was found transiently up-regulated in the people ascending to extremely high altitudes and was responsible for their adaptation [94]. The third SP1 target gene NR3C1 encodes for glucocorticoid receptor (GR) that controls in ammation and polymorphisms in the exon of NR3C1 are associated with the increased risk of HAPE in the Han Chinese population [94,95]. Lastly, the EGR1 regulates the SP1, an important in ammatory transcription factor. Recent studies have projected a negative correlation of EGR1 with the expression of interleukins expression IL-6 and IL-1 [96]. EGR1 has been known as a HAPE-induced gene in high altitude study [97]. Some studies also show that EGR1 remained elevated even during the recovery period from high altitude ascent. A recent study has identi ed an insertion/deletion (indel) polymorphism −742 indel in the promoter region of EPAS1 gene-speci c to Tibetans populations and not present in the lowlanders [98]. This insertion at −742 indel provided a new binding site for SP1 and was related to the activation of the EPAS1 promoter [91,98]. EPAS1 is well known for its contribution towards hypoxic response and is suggested to be responsible for the genetic adaptation of high-altitude hypoxia in Tibetans. The SP1 dependent transcriptional regulation is found to be responsible for cancer cell adaptation under hypoxic stress [91]. Recent pieces of evidence have shown that nuclear receptors are associated with metabolic adaptation in the Tibetans and genetic variants of nuclear receptors PPARG has been found highly associated with high-altitude hypoxia adaptation in native highlanders [70]. The PPARG is also under regulation of ESR1, which is an Estrogen nuclear receptor. The ESR1 protein is known to contribute to the activation of NOS3 and to endothelial production of nitric oxide and nitric oxide is a vasodilator that is found in high concentration in the high altitude native populations like Tibetans [101][102][103]. Estrogen receptors also help to accelerate the resolution of in ammation in M1 macrophage cells [104]. Hsa-miR-210-3p master regulate three overrepresented FFLs hsa-miR-210-3p-VEGFA-VHL, hsa-miR-210-3p-HIF1A-BNIP3 and hsa-miR-210-3p-HIF1A-LDHA. Both Hypoxia-inducible factor (HIF) and von Hippel-Lindau tumor suppressor protein (VHL) are hypoxia sensors that are found to control cellular responses to high altitude hypoxia [105]. Consortium 8 miRNAs and regulatory modules have been found to be associated with adipokines expression controlling metain ammation and probably regulating adaptation process at high altitude.
Examining the expression of both M1 and M2 cytokines in rat model under hypobaric hypoxic stress shows downregulation of M1 proin ammatory cytokine IL-6, while anti-in ammatory cytokines (TNFα, IL-1β) maintains expression close to the normoxia control after 72 hours. This shows that during high altitude hypoxic stress there is a shift in the M2 macrophages from M1 macrophages present in adipose tissue controlling in ammation state and providing protection from chronic in ammation. It also promotes anti-in ammatory response, tissue repair and insulin sensitivity.
Recent pieces of evidence have shown that orally transferred exogenous miRNAs derived from dietary substances could be absorbed through the intestine and remain active to regulate gene expression using the host machinery [24][25][26] . Hence a new concept of food-derived mRNAs appear as numerous reports of successful intervention of plant and animal source exogenous miRNAs have been utilized for augmented therapy in the mammalian system, demonstrating a new possible role of food beyond nutrition. E.g. plant miRNA miR2911 from honeysuckle (Lonicera periclymenum) against the in uenza-A virus in the mammalian system [106]. Viral miR-132 administration against Huntington's brain disease in a mammalian system [107] Plant miRNA miR-159a mimic sourced from Arabidopsis thaliana and soybean have been used to treat breast cancer [108]. hours exposure as compared control group but external miR-210-3p administration in high altitude hypoxia exposed mammalian have already found bene cial for protection against metabolic disorder like obesity and high altitude hypoxia adaptation. Detailed examination of the high-altitude studies also shows that 6 of the nine key miRNAs are relatively under expressed in low landers as compared to the native population. Contrastingly, hsa-miR-17-5p and hsa-miR-5582-5P are upregulated in non-acclimatize, low lander and healthy individuals. Hence hsa-miR-335-5p, hsa-miR-210-3p, hsa-miR-16-5p and hsa-miR-130a-3p have been considered for further studies. The uptake and bioavailability of miRNA are selective in nature and detailed examination of dietary miRNA studies has led to the identi cation of properties for screening potential dietary miRNA candidates. These properties are sequence conservation between interspecies miRNAs, high stability in body uids, and extracellular vesicle transport [65,109]. To know the miRNAs sequence conservation a phylogenetic study was conducted on the precursor sequences of differentially expressed miRNAs present in the high-altitude human population with high altitude mammals and non-mammals species. The 102 common miRNAs from Sus scrofa , Bos mutus, Bos taurus and non-mammals Gallus gallus were studied for phylogenetic relationship with human differentially expressed high altitude miRNAs hsa-miR-335-5p, hsa-miR-210-3p, hsa-miR-16-5p, hsa-miR-130a-3p that shares (>95%) of sequence similarity with the bta-miR-335, bta-mir-210, bta-mir-16-1, btamir-130a of Bos tauras species Table 1. Hence these miRNAs have the potential to be supplemented as dietary miRNAs to compensate for the downregulated complementary human miRNAs. To predict crosskingdom serum stability, sequence, and structural properties of these exogenous miRNAs have been checked considering the miR-168 as the gold standard [65]. Firstly, MFEI of all exogenous miRNAs is ≥ 0.85. Secondly, should ≤ 1 bulge in the sequences of these miRNAs and only 1-maximum 4 bulges are allowed. Finally, the sequences that <9 mismatches and ≤3 loops in the miRNAs stem-loop secondary structure Table 2. The results show that these miRNAs have the potential to be cross-kingdom stable in the serum of the mammalian system. There are three different types of miRNAs transfer and cellular uptake mechanism in a mammalian system and exogenous use of mammalian AGO2 protein to form miRNA-RISC complex and perform their functions similar to those of mammalian miRNAs is a common known mechanism [110,111]. Hence, docking studies were performed to predict 4 exogenous Bos tauras miRNAs and AGO2 binding e ciency. The result shows that all exogenous miRNAs have negative ACE, which means all binding reactions are spontaneous reactions but bta-miR-210 has the best ACE that is almost double the positive control (miR-168). The negative control is chosen as a randomly selected mammalian mir-20a of Cricetulus griseus (Chinese hamster). This shows that all 6 miRNAs bta-mir-16-1, bta-mir-130a, bta-miR-335, bta-mir-210 have the potential of dietary miRNA that can be used for early adaptation at high altitude. The complementary human miRNAs that are to be substituted are not only downregulated during high-altitude stress but some of them are found to regulate important functions during stress. E.g. miR-16 was found to be downregulated during HAPE and plays a key role in the disruption of the ion channels. Whereas, miR-26a has been found to inhibit the TGF-B signaling. Both regulations lead to the loss of cellular integrity and increased pulmonary pressure [112,113]. Lastly, miR-210 is termed as a master regulator of hypoxia and very recent evidence has shown that adipose-derived exosome miR-210 inhibits adipose browning to curb obesity and maintain body weight homeostasis [83].
The results also show that bta-miR-210 has the highest AGO2 binding e ciency and miR-210 chemically synthesized (agomiR) has proved to cure high-altitude obesity in the mammalian system [83]. This proves bta-miR-210 to be the most probable candidate for dietary miRNA therapeutics and can be provided with other 3 miRNA; bta-mir-16-1, bta-mir-130a, bta-miR-335 to compensate for complementary miRNAs. These miRNAs and their associated co-regulated module of TFs and target genes would help in better adaptation during rapid induction at high altitudes.