text:
Quantitative Prediction of Steatosis in Patients with Non-Alcoholic Fatty Liver by Means of Hepatic MicroRNAs Present in Serum and Correlating with Hepatic Fat

Non-alcoholic fatty liver disease (NAFLD) is the most prevalent form of chronic liver disease worldwide, but a reliable non-invasive method to quantify liver steatosis in primary healthcare is not available. Circulating microRNAs have been proposed as biomarkers of severe/advanced NAFLD (steatohepatitis and fibrosis). However, the use of circulating miRNAs to quantitatively assess the % of liver fat in suspected NAFLD patients has not been investigated. We performed global miRNA sequencing in two sets of samples: human livers from organ donors (n = 20), and human sera from biopsy-proven NAFLD patients (n = 23), both with a wide range of steatosis quantified in their liver biopsies. Partial least squares (PLS) regression combined with recursive feature elimination (RFE) was used to select miRNAs associated with steatosis. Moreover, regression models with only 2 or 3 miRNAs, with high biological relevance, were built. Comprehensive microRNA sequencing of liver and serum samples resulted in two sets of abundantly expressed miRNAs (418 in liver and 351 in serum). Pearson correlation analyses indicated that 18% of miRNAs in liver and 14.5% in serum were significantly associated with the amount of liver fat. PLS-RFE models demonstrated that 50 was the number of miRNAs providing the lowest error in both liver and serum models predicting steatosis. Comparison of the two miRNA subsets showed 19 coincident miRNAs that were ranked according to biological significance (guide/passenger strand, relative abundance in liver and serum, number of predicted lipid metabolism target genes, correlation significance, etc.). Among them, miR-10a-5p, miR-98-5p, miR-19a-3p, miR-30e-5p, miR-32-5p and miR-145-5p showed the highest biological relevance. PLS regression models with serum levels of 2–3 of these miRNAs predicted the % of liver fat with errors <5%.

Introduction
Non-alcoholic fatty liver disease (NAFLD), with a global prevalence of~25% [1,2], is one of the most important causes of liver disease [3] and will likely emerge as the leading cause of end-stage liver disease and transplantation soon. Its growing incidence in western countries is putting increasing pressure on healthcare systems.
The pathogenesis of NAFLD is considered multifactorial, and multiple mechanisms may be triggered during the course of the disease in each patient [4][5][6]. NAFLD involves, first, the accumulation of triglycerides (TG) in the liver, which may trigger a liver inflammatory response and damage. The histologic features range from isolated hepatic steatosis foci (non-alcoholic fatty liver; NAFL) to lobular inflammation and hepatocyte ballooning, which are hallmarks of non-alcoholic steatohepatitis (NASH). NASH very frequently precedes progressive fibrosis as a hallway to cirrhosis, liver-related and cardiovascular diseases, and hepatocellular carcinoma [4][5][6].
The current gold standard in NAFLD diagnosis and prognosis is the liver biopsy, where steatosis, lobular inflammation, ballooning, and fibrosis are semi-quantitatively scored. However, this is an invasive and not risk-free procedure, with frequent sampling errors (i.e., NAFLD does not affect the liver uniformly) and complications (pain, serious bleeding, injury to other organs, and, in rare cases, fatal outcomes). Moreover, histology interpretation is strongly influenced by pathologist training and perception [7]. However, while a number of non-or minimally invasive techniques (i.e., medical imaging) can rule out advanced fibrosis or cirrhosis, no completely reliable procedure to identify and score steatosis in early NAFLD patients is available [8].
Ultrasound has been used as a screening tool for fatty liver since it is non-invasive, inexpensive, and widely accessible in the clinical setting [9,10]. Nevertheless, the sensitivity and specificity of ultrasound imaging in accurately detecting the steatosis grade is controversial. Previous studies have shown wide ranges of sensitivity of 60-94% and specificity of 66-95% in detecting fatty liver. Other studies demonstrated that ultrasound is unreliable in detecting low-grade alterations in hepatic fat content. Moreover, ultrasound cannot reliably distinguish between fibrosis and steatosis. Finally, ultrasound does not provide reproducible quantitative information, and it can be influenced by the subjectivity of the examiner [9,10].
Other more reliable, non-invasive, imaging methods have other disadvantages that limit their use as a widely applicable screening tool for this prevalent disease. Computed tomography offers a semi-quantitative method for hepatic steatosis assessment, but it is distorted by an iron overload and involves exposure to radiation, which limits its use in longitudinal studies and in children. Imaging via magnetic resonance offers high accuracy for liver fat quantification without radiation exposure, but these scanners, in contrast to ultrasonographs, are costly, resource-demanding and frequently not applicable for certain patients (young children, claustrophobia, implanted electronic devices, metal implants, etc.) [9].
Consequently, there is an unmet need for reliable, accurate, simple, cost-affordable, quantitative, and non-or minimally invasive predictive NAFLD biomarkers in the clinical setting. Some biomolecules could meet this demand. MicroRNAs (miRNAs) are small, endogenous, noncoding RNAs involved in the control of many cellular pathways through the posttranscriptional regulation of target mRNAs, which are frequently released out of the cell into circulation. Liver miRNAs are definitively involved in the pathogenesis of NAFLD/NASH at various stages of the disease [11]. Consequently, several studies aimed at identifying circulating biomarkers for NAFLD have focused on miRNAs. Surprisingly, the miRNAs found in these independent studies show a very limited overlap, and only the liver specific miR-122 was unanimously claimed as a potential circulating NAFLD/NASH biomarker. Such a lack of agreement is likely to come from the bias in miRNA pre-selection, the reduced number of miRNAs examined, different technical approaches to quantify them, small study-size populations, and no biopsy-proven NAFLD diagnosis [12]. More recently, studies using large cohorts of patients and comprehensive miRNA profiling have overcome some of the previous shortcomings [13]. Nevertheless, relevant remaining questions are whether the serum miRNAs claimed as biomarkers (solely based on their biostatistical significance) do have a significant expression in the liver and show a quantitatively strong association with the steatosis grade, which is the early hallmark of NAFLD. It also needs to be determined whether the liver contributes significantly to the serum concentration of these postulated miRNAs and whether they have a potential implication in lipid metabolism and in the pathogenesis of NAFLD. Thus, it is conceivable that for an miRNA to be a reliable biomarker of hepatosteatosis, it should be liver-enriched, regulate lipid metabolism genes and correlate with the level of intrahepatic TGs. However, to our knowledge, the many studies investigating miRNAs as biomarkers of NAFLD have generally not addressed these logical premises.
Hepatosteatosis (lipid droplets in more than 5% of hepatocytes) is considered the primary event in the onset of NAFLD and is preceded by deregulated expression of key genes and proteins. Therefore, hepatosteatosis could also be associated with changes in miRNAs that regulate the expression of genes involved in lipid metabolism that, if released properly into the blood, could become early quantitative circulating NAFLD biomarkers.
Hence, our objective was, first, to comprehensively analyze the miRNAs expressed in steatotic human livers and examine their correlation with the concentration of liver TGs; and second, to analyze all human serum miRNAs present in NAFLD patients histologically diagnosed with different steatosis grades and examine their correlation with the % of liver fat in the stained liver biopsies-all of this to identify common liver and serum deregulated miRNAs in steatotic patients, score their relevance, and develop a serum-based model with the most relevant miRNAs able to predict the % of liver fat in patients with early NAFLD.

Human Liver and Serum MicroRNA Sequencing and Building of Steatosis Predictive Models
Comprehensive microRNA sequencing of 20 liver samples from cadaveric organ donors resulted in a list of 418 miRNAs detected in at least 18 samples (90%). From these 418 miRNAs, 89 (18%) showed a significant Pearson correlation with the concentration of liver TGs. Principal component analysis (PCA) showed an association between the scores of the first principal component (PC), explaining 31.4% of the variance, and the TG content ( Figure 1A).
(solely based on their biostatistical significance) do have a significant expression in the liver and show a quantitatively strong association with the steatosis grade, which is the early hallmark of NAFLD. It also needs to be determined whether the liver contributes significantly to the serum concentration of these postulated miRNAs and whether they have a potential implication in lipid metabolism and in the pathogenesis of NAFLD. Thus, it is conceivable that for an miRNA to be a reliable biomarker of hepatosteatosis, it should be liver-enriched, regulate lipid metabolism genes and correlate with the level of intrahepatic TGs. However, to our knowledge, the many studies investigating miRNAs as biomarkers of NAFLD have generally not addressed these logical premises.
Hepatosteatosis (lipid droplets in more than 5% of hepatocytes) is considered the primary event in the onset of NAFLD and is preceded by deregulated expression of key genes and proteins. Therefore, hepatosteatosis could also be associated with changes in miRNAs that regulate the expression of genes involved in lipid metabolism that, if released properly into the blood, could become early quantitative circulating NAFLD biomarkers.
Hence, our objective was, first, to comprehensively analyze the miRNAs expressed in steatotic human livers and examine their correlation with the concentration of liver TGs; and second, to analyze all human serum miRNAs present in NAFLD patients histologically diagnosed with different steatosis grades and examine their correlation with the % of liver fat in the stained liver biopsies-all of this to identify common liver and serum deregulated miRNAs in steatotic patients, score their relevance, and develop a serumbased model with the most relevant miRNAs able to predict the % of liver fat in patients with early NAFLD.

Human Liver and Serum MicroRNA Sequencing and Building of Steatosis Predictive Models
Comprehensive microRNA sequencing of 20 liver samples from cadaveric organ donors resulted in a list of 418 miRNAs detected in at least 18 samples (90%). From these 418 miRNAs, 89 (18%) showed a significant Pearson correlation with the concentration of liver TGs. Principal component analysis (PCA) showed an association between the scores of the first principal component (PC), explaining 31.4% of the variance, and the TG content ( Figure 1A).  A PLS model was built, providing a statistically significant RMSECV = 421 (permutation test p-value = 0.01). The recursive feature elimination (RFE) algorithm was employed to find the subset miRNA providing the lowest RMSECV in predicting the %TG. From results depicted in Figure 2A, a subset of 50 miRNAs was selected, providing an RMSECV = 240 (57% of the error obtained using the whole miRNA set).
A PLS model was built, providing a statistically significant RMSECV = 421 (permutation test p-value = 0.01). The recursive feature elimination (RFE) algorithm was employed to find the subset miRNA providing the lowest RMSECV in predicting the %TG. From results depicted in Figure 2A, a subset of 50 miRNAs was selected, providing an RMSECV = 240 (57% of the error obtained using the whole miRNA set). The same strategy was followed for the analysis of serum miRNA profiles. Whole miRNome sequencing of 23 serum samples from NAFLD patients resulted in a list of 351 miRNAs that were detected in at least 19 serum samples (82.5%). From these 351 miRNAs, 51 (14.5%) showed a significant Pearson correlation with the % of fat in the paired liver biopsies of these patients. PCA showed an association between the scores of the first PC, explaining 22.0% of the variance, and the % of fat ( Figure 1B). A PLS model was built using the complete set of retained miRNAs, providing a statistically significant RMSECV = 6 (permutation test p-value = 0.03). Results by RFE enabled the identification of a subset of 50 miRNAs, providing an RMSECV = 3 (50% of the error initially obtained using the whole miRNA set) ( Figure 2B). The comparison of the two subsets of miRNAs selected from the liver and serum datasets showed that 19 miRNAs (listed in Table 1) were coincident.  The same strategy was followed for the analysis of serum miRNA profiles. Whole miRNome sequencing of 23 serum samples from NAFLD patients resulted in a list of 351 miRNAs that were detected in at least 19 serum samples (82.5%). From these 351 miRNAs, 51 (14.5%) showed a significant Pearson correlation with the % of fat in the paired liver biopsies of these patients. PCA showed an association between the scores of the first PC, explaining 22.0% of the variance, and the % of fat ( Figure 1B). A PLS model was built using the complete set of retained miRNAs, providing a statistically significant RMSECV = 6 (permutation test p-value = 0.03). Results by RFE enabled the identification of a subset of 50 miRNAs, providing an RMSECV = 3 (50% of the error initially obtained using the whole miRNA set) ( Figure 2B). The comparison of the two subsets of miRNAs selected from the liver and serum datasets showed that 19 miRNAs (listed in Table 1) were coincident.

Biological Significance of the 19 miRNAs Identified Both in Liver and Serum That Show Correlation with Liver Steatosis
First, most of the selected miRNAs matched the guide strand of their miRNA duplex ( Table 1). The guide strand is retained in Ago proteins and stably forms the RISC complex. The other strand, known as the passenger strand (or miRNA *), is discarded and usually degraded. Nevertheless, in some cases, both arms of the duplex give rise to functional mature miRNAs that are loaded into Ago proteins. As our goal was to identify circulating miRNA biomarkers, the more abundant guided strands were preferred. Secondly, we assessed the trend and strength of the association of each miRNA with liver TGs and % of fat. We unexpectedly found that for most of the selected miRNAs (80%), the positive or negative association with steatosis in the liver usually translated into an association with the opposite trend in serum (Table 1, Figure 3). Thus, most miRNA that increased with TG levels in the liver, steadily decreased in serum with the % of liver fat, and vice versa. Results suggest that, for most miRNAs, hepatosteatosis influences the dynamics of miRNA release rather than miRNA biogenesis.
Note: *: Ratio of mean values in liver and serum samples.

Biological Significance of the 19 miRNAs Identified Both in Liver and Serum That Show Correlation with Liver Steatosis.
First, most of the selected miRNAs matched the guide strand of their miRNA duplex ( Table 1). The guide strand is retained in Ago proteins and stably forms the RISC complex. The other strand, known as the passenger strand (or miRNA *), is discarded and usually degraded. Nevertheless, in some cases, both arms of the duplex give rise to functional mature miRNAs that are loaded into Ago proteins. As our goal was to identify circulating miRNA biomarkers, the more abundant guided strands were preferred. Secondly, we assessed the trend and strength of the association of each miRNA with liver TGs and % of fat. We unexpectedly found that for most of the selected miRNAs (80%), the positive or negative association with steatosis in the liver usually translated into an association with the opposite trend in serum (Table 1, Figure 3). Thus, most miRNA that increased with TG levels in the liver, steadily decreased in serum with the % of liver fat, and vice versa. Results suggest that, for most miRNAs, hepatosteatosis influences the dynamics of miRNA release rather than miRNA biogenesis.  It was also important to consider the absolute value and significance of the Pearson correlation coefficient. Some serum miRNAs selected by the PLS-RFE model (i.e., 148a-3p, 660-5p, 335-5p and 17-3p) did not show a significant correlation (Table 1), and this could limit their convenience as individual predictive circulating biomarkers. To check the biological relevance of the selected miRNAs, we used the DIANA-microT-CDS [14] and the mirDIP algorithms [15] against the collection of lipid metabolism genes, as specified by REACTOME Id: R-HSA-556833, containing 749 proteins from 730 genes. Many of the selected miRNAs were predicted to regulate numerous lipid metabolism genes, but some of them (e.g., 191-5p, 769-5p, 660-5p, 30a-3p, and 136-3p) had a very limited impact on these genes (Table 1). Finally, we assessed the abundance of these miRNAs in the liver. It is conceivable that a circulating miRNA biomarker of hepatosteatosis will be significantly expressed in and released from the liver. Moreover, it is likely that liver-specific miRNAs will be diluted in the blood, resulting in high liver level/serum level ratios. The more tissues contributing to the serum concentration of a given miRNA, the lower the liver/serum ratio will be. Liver-enriched miRNAs 122-5p and miR-192-5p exhibit liver/serum ratios between 300 and 3000. We observed low liver/serum ratios with several of the selected miRNAs, such as let-7d-5p, let-7b-5p, 191-5p or 140-3p (Table 1).

PLS Model Based on Only 2-3 Serum miRNAs to Predict % of Fat in NAFLD Patients
A series of PLS models including two or three miRNAs from the subset of 19 common miRNAs (quantified in serum samples) were investigated by exhaustive feature selection. The performance for the accurate prediction of liver fat % was again evaluated. The combinations of 2-3 miRNAs with higher biological significance and more accurate prediction of fat % (lower RMSECV) were finally selected (Table 2, Figure 4). Therefore, according to these models, it would be possible to predict the % of fat in the liver of NAFLD patients by quantifying 2-3 miRNAs in the serum with an error of around 4-5%. This could become a very useful tool to screen early NAFLD in the general population and stratify patients according to their liver fat %.

Discussion
A group of experts has recently proposed a change in nomenclature from non-alcoholic fatty liver disease (NAFLD) to metabolically associated fatty liver disease (MAFLD). The new term better reflects the pathophysiology of NAFLD as a metabolically driven disease, shifting from "exclusion" criteria to "positive" diagnosis criteria. One important additional difference is that NAFLD diagnosis is based primarily on histological criteria, being steatosis-graded by anatomical pathologists. However, for MAFLD, steatosis may be diagnosed non-invasively with either imaging techniques or biomarkers [16]. In this regard, ultrasonography, contrary to other imaging techniques, is widely available and

Discussion
A group of experts has recently proposed a change in nomenclature from non-alcoholic fatty liver disease (NAFLD) to metabolically associated fatty liver disease (MAFLD). The new term better reflects the pathophysiology of NAFLD as a metabolically driven disease, shifting from "exclusion" criteria to "positive" diagnosis criteria. One important additional difference is that NAFLD diagnosis is based primarily on histological criteria, being steatosis-graded by anatomical pathologists. However, for MAFLD, steatosis may be diagnosed non-invasively with either imaging techniques or biomarkers [16]. In this regard, ultrasonography, contrary to other imaging techniques, is widely available and could be the first-line diagnostic tool for assessing suspected hepatic steatosis or a screening tool for asymptomatic NAFLD [17]. However, while ultrasonography could be reliable at detecting moderate steatosis, its sensitivity is thought to be poor when <20-30% of hepatocytes are steatotic [17]. For instance, in a comparative study of different imaging techniques, ultrasonography showed the lowest sensitivity (65%) and specificity (77%) in detecting ≥5% histologically defined hepatic steatosis [18]. In addition, ultrasonography usually reports steatosis on a semiquantitative scale (normal, mild, moderate, and severe). Thus, the complementation of widely available imaging techniques along with circulating biomarkers could increase accuracy in steatosis diagnosis, and with a quantitative score.
Several classic biomarkers are included in algorithms intended to score steatosis. However, most of these biomarkers are nonspecific enzymes unrelated to liver TG accumulation. The hepatic steatosis index, NAFLD liver fat score, and fatty liver index include typical liver enzymes such as hepatic transaminases or γGT. However, it is well known that normal enzyme levels can be found in NAFLD patients representing the entire spectrum of disease severity, from simple steatosis to advanced NASH [19].
The search for novel useful biomarkers for early NAFLD that could routinely and economically be applied to patients to diagnose, screen, and/or monitor the extent and progress of NAFLD led us to examine in detail biomolecules that originated in the liver, are present in patient's sera, and could account for the % of fat in the liver. By using two separate cohorts, we identified, from the large number of miRNAs present in liver tissue and sera of individuals with a variable grade of steatosis, those having a stronger association with the % of liver fat and mechanistically involved in lipid metabolism regulation.
A key feature making serum miRNAs excellent potential biomarkers is their stability. They are resistant to different external insults because they do not circulate as free RNA; instead, they are encapsulated in membranous vesicles (micro-vesicles, exosomes, apoptotic bodies), complexed to RNA-binding proteins (e.g., Ago2) or associated with lipoproteins, which protect them from endogenous RNases. Another aspect making miRNA molecules attractive biomarkers is that the RT-qPCR technique used for their detection is extremely sensitive and cost-effective [20].
Several previous studies have also investigated miRNA as biomarkers for NAFLD. In a previous study, we attempted to validate many of these miRNAs postulated in previous studies, but the rate of success was very low. Our results suggested that measurement design, technical and methodological differences, the pre-selection of miRNA to be measured, and the reduced number of miRNAs examined likely limited reproducibility [12]. In the present study, a comprehensive miRNAseq encompassing all of the expressed miRNAs in serum and liver tissue was performed. Another important bias is that in some studies, NAFLD diagnosis is not biopsy-proven, thus opening the possibility that false positive and negative patients could be inadvertently included. Our study is the first including not only biopsy-proven patients but also patients with an accurate quantification of steatosis either by biochemical methods (µg TG/mg liver protein) or by image quantifications (based on the % of the area occupied by lipid droplets in the hematoxylin-eosin stained liver biopsy). As an important contribution to minimize subjective biasing in the score of steatosis, we had high-quality and sufficiently large biopsies analyzed by the same pathologist in a blind manner; that is, she was not aware of the clinical features and diagnosis of the patients at the time of examination. On top of that, steatosis was estimated not only visually but also by computer analysis of images taken from various fields of hematoxylin-eosin stained tissue sections. This image quantification was based on previous studies [21] and validated before being applied in this work by comparing computer data analysis and the lipid content (biochemically determined) from the same liver samples.
Out of the large number of miRNAs identified in liver tissue and serum, a much smaller number fulfilled the criteria we set ourselves, i.e., (a) be liver-enriched, (b) be involved in the regulation of lipid metabolism genes, and (c) correlate significantly with the level of intrahepatic lipids. Only a few of the model-selected miRNAs met these criteria.
The relevance of some of the miRNAs selected after our regression modeling and biological fine-tuning was reinforced by previous studies searching for miRNA biomarkers. Thus, miR-23a-3p and miR-19a-3p were found to be deregulated (with an opposite trend) in the serum of patients with NASH [22]. Similarly, miR-19a and b were both upregulated in the serum of simple steatosis NAFLD patients [23], whereas in a different study, miR-30a-3p was also substantially upregulated in NAFLD patients, and its suppression attenuated hepatic steatosis in HepG2 cells [24]. In an animal model, miR-145-5p increased in steatotic liver in response to a high-fat diet [25]. Moreover, ingestion of a high-fat high-saturated meal was able to deregulate miR-145-5p in human serum [26].
Several of the miRNAs proposed as steatosis biomarkers have demonstrated links to key important lipid metabolism pathways. Thus, miR-98-5p inhibits PGC-1β mRNA expression in the liver, and PGC-1β plays an important role in the regulation of hepatic lipogenesis by activating genes such as Scd1 and Fas [28]. The decrease in liver miR-98-5p observed in our study would likely result in more PGC-1β and lipogenesis.
In a different study, miR-32-5p targeted KLF3, which in turn represses lipogenic regulators such as PPARγ and C/EBPα [29]. In this case, the decrease in liver miR-32-5p observed in our study will result in higher KLF3 and lower lipogenesis as a likely compensatory response.
Another recent study showed that hepatocytes release let-7b-5p after stimulation with palmitic acid by a TGFβ-dependent mechanism. Moreover, let-7b-5p overexpression increased hepatocyte fatty acid transport [30]. In our study, this miRNA was upregulated both in the liver and serum.
Nevertheless, the potential utility of the proposed novel miRNA biomarkers might be limited by other factors, such as the complex contribution of other co-morbid pathologies, particularly obesity, T2D and dyslipidemia. Further improvements to give confidence to our findings would include performing a paired analysis of liver and serum miRNAs from the same cohort of patients, as well as running a prospective study to monitor the progression of the disease by miRNA analysis, medical imaging, and whenever possible, liver biopsy. This clinical study is currently ongoing.
Moreover, due to the limited sample size of the two study cohorts, further efforts will be required to validate our results in much larger sample size cohorts. However, if these results are confirmed in future studies, patients with NAFLD may benefit from diagnostic or surveillance programs based on minimally invasive miRNA biomarkers.

Patients and Human Samples
Twenty human liver samples were obtained from organ donors who were primarily assigned to liver transplantation, but their livers failed to be transplanted and were donated to research (Biobank Hospital La Fe, Valencia, Spain). Histological analyses were not performed on these livers, but we did perform biochemical analyses of TG and total liver lipids and found a wide range of variation (Table 3). Anthropometric and analytical characteristics of these liver donors have been described elsewhere [31] and are summarized in Table 3. Cadaveric organ donors comprised 14 men and 6 women, and their causes of death included cerebrovascular accident (either hemorrhagic or ischemic, n = 15) or cranioencephalic traumatism (n = 5). All of them were not harboring any infectious disease and tested negative for human immunodeficiency and hepatitis viruses. Livers were obtained in the operating room and transported in a cold preservation solution (University of Wisconsin solution). After the reception, and once rejected for transplantation, livers were immediately cut into small pieces, frozen in liquid nitrogen, and stored at −80 • C until use. In this cohort of donors, the hepatic levels of intrahepatic TG and total lipids covered a wide range (10-20-fold variation). The predominant lipids in the liver were TG and consequently, the correlation between total lipids and TG was very high (R Pearson = 0.91). We also observed significant variability in serum bilirubin and transaminases (AST and ALT) (Table 3). However, they did not show any correlation with liver TG levels, thus ruling them out as collinear confounding variables. It is important to remark that biochemical data from these donors were not necessarily obtained in standard conditions, and some parameters might be influenced by other factors, such as the timing of blood extraction and cerebral death.
Human sera were collected from 23 well-characterized NAFLD patients, covering all grades of steatosis as confirmed by percutaneous liver biopsy. The NAFLD patients did not consume alcohol regularly (less than 20 g/day), and other potential causes of liver disease (viral, autoimmune, hepatotoxic drugs, iron overload, Wilson's disease, etc.) were excluded. Hematoxylin-eosin and Masson's trichrome-stained paraffin-embedded liver biopsy sections were examined and interpreted by the same experienced hepatopathologist (J.P.-R.), who was unaware of the patients' clinical data. All liver biopsies examined showed more than 10 complete portal tracts. Steatosis along with ballooning, lobular inflammation and fibrosis were assessed as outlined by Brunt et al. [32]. Disease severity was scored according to both NAS [32,33] and SAF [34,35] systems. Steatosis, as % of fat in the biopsy, was automatically quantified by image analysis at 20× resolution of hematoxylin-eosin stained sections in MATLAB 2021b (Mathworks Inc., Natick, MA, USA) using in-house written scripts and the Image Processing Toolbox (MATLAB), following a strategy similar to that previously described by Munsterman et al. [21]. Anthropometric and analytical characteristics of these NAFLD patients have been described elsewhere [12] and are summarized in Table 4. Anticoagulated venous blood was extracted between 8 and 10 am, after overnight fasting, by the time of liver biopsy. Blood was collected in siliconized tubes and centrifuged at 2500× g for 10 min. Serum samples were stored at −80 • C.

Quantification of Intrahepatic TG and Total Lipids
Homogenates from human livers were extracted with a methanol-chloroform mixture and evaporated under nitrogen as described elsewhere [36]. TG and total lipid concentrations were measured in the lipid residue using colorimetric kits from Spinreact (Gerona, Spain) based on the GPO-POD enzymatic method (#1001311) and the sulfo-phospo-vanillin reaction (#1001270), respectively. Protein concentration was determined using the Protein Assay Kit from Bio-Rad Laboratories (Madrid, Spain).

RNA Isolation, Small RNA Library Preparation and miRNAseq Analysis
Total RNA was isolated from liver and serum samples using the miRNeasy mini kit (Qiagen Westburg BV, Leusden, The Netherlands) according to the manufacturer's protocol, followed by a DNAse I treatment (Qiagen Inc., Venlo, The Netherlands). RNA concentration and quality were measured using a BioAnalyzer system (Agilent Technologies, Breda, the Netherlands). Starting from total RNA samples, small RNAs were size selected and ligated for sequencing following the TruSeq Small Prep Kit Preparation (15004197 Rev. D, Illumina, Eindhoven, The Netherlands). Samples were sequenced on the HiSeq 2000 (Illumina) in single-end 50 bp. Sequencing data have been deposited on ENA (European Nucleotide Archives) under the accession number PRJEB53387. Small RNA reads were first trimmed from the 5 adapter sequence, and only the post-trimmed reads of 16 to 35 bp were selected. All trimmed reads were then mapped using Patman [37] against the human precursor database from miRbase (release 11), allowing no mismatches or gaps. Patman outputs were parsed to obtain all read mapping to each miRNA arm (5 and 3 ). Normalization and expression analysis were performed using the R package DESeq2 [38].

Bioinformatics Analysis and Modeling
The analysis initially included miRNAs commonly detected in both liver and serum samples. A data clean-up step was carried out to improve the robustness of the analysis. Accordingly, miRNAs detected in <90% or <82.5% of the liver and serum samples, respectively, as well as those miRNAs showing median reads higher in serum than in liver, were removed, leaving a set of 154 miRNAs for further analysis.
Two independent multivariate partial least squares (PLS) regression methods were built to predict: (a) the concentration of TG in human liver based on liver miRNAs and (b) the steatosis level (% of fat in biopsy images) in human liver based on serum miRNAs. PLS was selected for multivariate regression because of the expected multicollinearity among the covariates (miRNAs) and the large number of miRNAs (i.e., 154) compared to the smaller sample size (i.e., [20][21][22][23]. The selection of the number of PLS latent variables was based on results obtained by leave-one-out cross validation (LOOCV). The assessment of the significance of the figures of merit (e.g., root mean square error of cross validation, RMSECV), was carried out by permutation testing (500 permutations), where the p-value was estimated as the fraction of the permuted models, providing a better (i.e., lower RMSECV) estimate than that estimated using the original y vector.
The recursive feature elimination (RFE) algorithm was employed to find the miRNA most relevant in predicting the target variable (TG concentration or % of fat assessed in the biopsy images). Here, a PLS full model was created, and a measure of variable importance (VIP score) was computed to rank the predictors from most important to least. At each stage of the search, the least important predictor was eliminated prior to rebuilding the model. At each iteration, the PLS model performance was assessed using LOOCV, and a new RMSECV was reported. The analysis of the evolution of the obtained RMSECV during the RFE searches was used for the identification of two subsets of miRNAs for each model providing the optimal model performance. The analyses of the two datasets by PLS-RFE enabled the identification of 50 miRNAs in each of the datasets associated with the liver TG concentration and biopsy fat%, and 19 among them were commonly selected in both cases.
Finally, PLS regression models to predict steatosis grade with 2 or 3 miRNAs were built. The low number of miRNAs included in this step allowed us to evaluate all possible combinations of the list of 19 miRNAs, taken 2 or 3 at a time (i.e., 171 and 969, respectively). The evaluation of these models was based on the RMSECV and on the biological relevance of selected miRNAs in lipid metabolism and NAFLD-associated pathways assessed through an in silico search of mRNA targets.

Conclusions
In this proof-of-concept study, we demonstrated that it is possible to identify serum miRNAs, with biological significance, that correlate with the % of liver fat in NAFLD patients. Moreover, our findings allowed us to suggest that a panel of biomarker miRNAs, identified according to the strategy of this study, could become a promising tool to detect and stratify patients with NAFLD through a minimally invasive procedure. Thus, algorithms based on circulating miRNAs could facilitate a "next level" quantitative diagnosis and clinical follow-up of NAFLD patients, even at early stages of the disease.


Statistical and functional convergence of common and rare genetic influences on autism at chromosome 16p

The canonical paradigm for converting genetic association to mechanism involves iteratively mapping individual associations to the proximal genes through which they act. In contrast, in the present study we demonstrate the feasibility of extracting biological insights from a very large region of the genome and leverage this strategy to study the genetic influences on autism. Using a new statistical approach, we identified the 33-Mb p-arm of chromosome 16 (16p) as harboring the greatest excess of autism’s common polygenic influences. The region also includes the mechanistically cryptic and autism-associated 16p11.2 copy number variant. Analysis of RNA-sequencing data revealed that both the common polygenic influences within 16p and the 16p11.2 deletion were associated with decreased average gene expression across 16p. The transcriptional effects of the rare deletion and diffuse common variation were correlated at the level of individual genes and analysis of Hi-C data revealed patterns of chromatin contact that may explain this transcriptional convergence. These results reflect a new approach for extracting biological insight from genetic association data and suggest convergence of common and rare genetic influences on autism at 16p.

Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. further. Nature Genetics is committed to improving transparency in authorship. As part of our efforts in this direction, we are now requesting that all authors identified as 'corresponding author' on published papers create and link their Open Researcher and Contributor Identifier (ORCID) with their account on the Manuscript Tracking System (MTS), prior to acceptance. ORCID helps the scientific community achieve unambiguous attribution of all scholarly contributions. You can create and link your ORCID from the home page of the MTS by clicking on 'Modify my Springer Nature account'. For more information please visit please visit www.springernature.com/orcid.
We look forward to seeing the revised manuscript and thank you for the opportunity to review your work.

4
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Remarks to the Author:
Authors introduce a stratified-pTDT (S-pTDT), which estimates transmission in parent-child trios of PRS constructed from regions/blocks of the genome. Authors tested whether S-pTDT could identify any regions of the genome with transmission of ASD polygenic risk significantly over or under genome-wide expectation for large blocks of the genome. Authors estimated transmission two large trio samples and the transmission of regional polygenic risk for ASD is correlated between the two samples. Partitions with large S-pTDT z-scores cluster on 16p (0-33Mb).
The 16p does not contain a genome-wide significant locus for ASD, authors still sought to determine whether the S-pTDT signal at 16p could be explained by one or a small number of common variant associations -single driving locus in the region was not found and authors verified that the overtransmission of ASD polygenic risk at 16p is not driven by CNV carriers in their data.
16p is gene rich and many of the genes are only expressed in the brain. Gene density, density of brainspecific genes or density of constrained genes, based on authors' inquiries, cannot explain the region's degree of polygenic over-transmission.
Interestingly, authors observe across independent cohorts that increased 16p ASD PRS is associated with an average decrease in expression of brain expressed genes within the 16p region. Also, in vitro deletion of the 16p11.2 locus was associated with decreased average expression of 200 neuronally expressed genes on chromosome 16p. Furthermore, authors suggest that the 16p region may have increased within-region chromatin contact, which could explain the apparent non-independence of genetic and expression variation at mega-base scale. Authors hypothesize that this diffusely elevated within-region contact at 16p could facilitate the influence of regional polygenic effects on gene expression across 16p, via complex distal regulatory interactions. Lastly authors conclude that the 16p11.2 CNV has increased physical interaction with the telomeric region and the 3D conformation of 16p may mediate convergent ASD-related genetic effects on gene expression via regulatory interactions across mega-bases of separation. Based on these observations authors present the "Integrative model of ASD liability at 16p". Comments 1. The samples studied are large and impressive, the analyses are transparent, novel and sound and the model is interesting. The main weakness is how modest the mean expression effects are for the 16p region, particularly when compared to the decrease in gene expression associated with heterozygous gene deletion (16p11.2 CNV).
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. 2. The deletion described and studied in this manuscript is the 16p11.2 proximal deletion explaining up to 1% of autistic cases. It is particularly interesting that distal to this locus is another recurrent deletion, 16p11.2 distal deletion, not mentioned in the manuscript. That deletion, the 16p11.2 distal deletion confers high-risk of very similar phenotypes (including autism, cognitive impairments and obesity). If that deletion, the 16p11.2 distal deletion, also affects gene expression on 16p in similar manner as the 16p11.2 proximal deletion and the 16p ASD PRS, the story would be more convincing. Thus, my recommendation is to include as well data on the 16p11.2 distal deletion in the manuscript. In this study, Weiner et al. investigate the feasibility of extracting biological insight from a large genomic region and understanding how it is associated with risk for autism spectrum disorder. They identified the 33 Mb short arm of chromosome 16 as harboring the greatest excess of common polygenic risk for ASD. Analysis of bulk and single-cell RNA-sequencing data from post-mortem human brain samples revealed that common polygenic risk for ASD within 16p is associated with decreased average expression of genes throughout this 33-Mb region.
They subsequently use isogenic neuronal cell lines with CRISPR/Cas9-mediated deletion of 16p11.2 to show that the deletion is also associated with depressed average gene expression across the short arm of chromosome 16. The effects of the rare deletion and diffuse common variation were correlated at the level of individual genes. Their results also suggest that very dense 3D chromatin contact within the short arm of chromosome 16 may coordinate genetic and transcriptional disease liability across this region.

6
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. This study focuses on a large genomic segment rather than on specific genes. It is original and advances the field. The claims are supported by the data. In particular, the authors provide a rather convincing link between the transcriptional effects of rare and common variants implicated in ASD.
My major comment is that there is a lack of functional characterization of this large group of genes on the short arm of 16p. If these genes are co-regulated, this would imply that they are implicated in shared functional modules. The study would benefit from a functional characterization: i.e., Are genes within this genomic block enriched in well-known functional modules? Or in modules identified by contrasting gene expression in the brains of individuals with ASD and controls?
Specific comments: In the 2nd paragraph of the results: "... we constructed stratified PRS from adjacent blocks of SNPs, yielding 2,006 (often overlapping) partitions collectively covering the whole genome (median number of SNPs per block: 3,000, minimum length: 4.3Mb, maximum length: 52.9Mb, median length: 11.7Mb, Supplementary Figure 3, Methods)..." It is unclear how they defined these genomic blocks of very different sizes. Was it based on the number of LD blocks, the number of genes? Or completely random? Figure 1E. The number of blocks removed stops at n=25. Is that just because there was no more effect? One would expect that the SE of the pTDT would get larger as more blocks are removed, but that doesn't appear to be the case. It also seems like there is a trend that may become significantly protective at one point. What would happen beyond =25? In other words, once authors remove the most over transmitted blocks, are there protective blocks in the 16p11.2 short arm?
As a sensitivity analysis, the authors performed an analogous analysis using a cohort of ADHD trios and an external ADHD GWAS and they did not replicate the finding in ADHD. One could argue that ADHD may be the worst condition to perform such a sensitivity analysis since the PRS doesn't explain much variance. Schizophrenia would appear to be much more relevant. The PRS is more robust, cohorts are larger, and the 16p11.2 locus is associated with schizophrenia.
The relationship between gene density and over-transmission is an important point and should be represented in a figure in the main text.
The authors asked whether, on average, the 200 neuronally expressed genes on 16p were differentially Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. expressed in response to the 16p11.2 deletion. Genes on 16p had significantly lower expression in the deletion lines. The deletion's effect on 16p genes differed from the effect on all other 8,533 neuronally expressed genes in the genome (P = 0.02 -somewhat of a trend -), whose expression was not, on average, changed by the deletion (P = 0.43).
Can the authors provide information on the non-neuronally expressed genes on the short arm? Would these represent a better "control group", providing a stronger contrast? Shouldn't the Y-axis of figure 2C represent "fold change" instead of t-stat?
Increased ASD PRS within 16p was associated with decreased expression in glutamatergic neurons of genes through the 16p region. Do the authors observe an increase in the variance of gene expression? In other words, is this a simple shift in mean expression, or do results also suggest that there may also be some genes with an increase in expression?
Authors show that the CNV-telomeric contacts (n = 291 100kb x 100kb contacts) are 2.9x more frequent than contacts between distance-matched control regions on 16p (n = 1,808 100kb x 100kb contacts, P < 1e-10). However, the 16p11.2 region is 30MB away from the telomeric region. I don't see, therefore how they can test distance-matched regions on the short arm. Distance-matched regions would only be found downstream of the 16p11.2 locus on the long arm.
Authors suggest that the entire short arm may represent a group of co-regulated genes involved in ASD. It would be necessary to demonstrate that this is the case and test the enrichment of these genes in known functional modules (in health and disease. i.e.). For example, are these gene enriched differentially expressed genes obtained by contrasting gene expression in the brains of individuals with ASD and controls. This data is available, and the analysis should be straightforward.

Response to reviewers
We thank the reviewers for their thoughtful comments and have responded point-by-point below. We have also highlighted any corresponding changes in the main text.

Reviewer #1:
Comment #1.1: "The deletion described and studied in this manuscript is the 16p11.2 proximal deletion explaining up to 1% of autistic cases. It is particularly interesting that distal to this locus is another recurrent deletion, 16p11.2 distal deletion, not mentioned in the manuscript. That deletion, the 16p11.2 distal deletion confers high-risk of very similar phenotypes (including autism, cognitive impairments and obesity). If that deletion, the 16p11.2 distal deletion, also affects gene expression on 16p in similar manner as the 16p11.2 proximal deletion and the 16p ASD PRS, the story would be more convincing. Thus, my recommendation is to include as well data on the 16p11.2 distal deletion in the manuscript." We thank the reviewer for this thoughtful comment about the possibility that other ASD-associated CNVsespecially those located on 16p -may confer disease liability through similar mechanisms to the proximal 16p11.2 studied here. The reviewer notes the potential relevance of the 16p11.2 distal deletion, and we agree this proximate and disease-associated CNV would be interesting to investigate. However, as far as we can tell, there are no published whole-genome RNA-sequencing datasets of the 16p11.2 distal deletion, including in the provided reference of Sønderby et al. It is also materially infeasible for us to generate additional isogenic distal deletions at this time. Without RNA-sequencing, we are unable to evaluate our model for this deletion.
That said, we are very interested in whether other neuropsychiatric CNVs confer disease liability through similar mechanisms to those described for the proximal 16p11.2 CNV. We are actively testing this hypothesis in our research group, and we look forward to sharing our findings with the community in the future. Finally, we clarified in the manuscript that we analyzed the proximal and not the distal deletion (page 3).

Reviewer #3:
Comment #3.1: "My major comment is that there is a lack of functional characterization of this large group of genes on the short arm of 16p. If these genes are co-regulated, this would imply that they are implicated in shared functional modules. The study would benefit from a functional characterization: i.e., Are genes within this genomic block enriched in well-known functional modules? Or in modules identified by contrasting gene expression in the brains of individuals with ASD and controls?" We thank the reviewer for raising this important question -we are also extremely eager to understand the aggregated downstream functional consequence of genetic variation in the region.
We first performed gene ontology (GO) analysis to evaluate enrichment of genes on 16p in annotated biological pathways (http://geneontology.org/). We used the same 17,909 genes from the gene density analysis as reference genes. We tested for enrichment of all genes on 16p (midpoint < 32,000,000 bp, n = 432 genes) across three classes of annotations: biological process, molecular function, and cellular component.
The GO analysis for molecular function and cellular component returned multiple bonferroni-significant enrichments: multiple lipid/fatty acid pathways (Fatty-acyl-CoA synthase activity, Butyrate-CoA ligase activity, Medium-chain fatty acid-CoA ligase activity, >20x enrichment for each), and hemoglobin complex (19x enrichment). The lipid/fatty acid pathways return an enrichment because there are 5 acyl-CoA-synthase genes located within 500kb of each other on 16p around Mb 20. Similarly, the hemoglobin complex pathway returns an enrichment because 4 hemoglobin subunits are clustered together within 100kb of each other at the start of chromosome 16. These examples raise a critical point: since functionally similar genes are often clustered together in the genome (Andrews et al. 2015 Genome Research), a gene set enrichment signal will be dominated by whichever functional cluster of genes happens to be located within the region of interest. Thus, we do not believe that canonical gene set enrichment approaches are suited to regional enrichment analysis. That said, it is also possible that decreased expression across 16p does not exert direct phenotypic effect, but instead propagates to interact with gene/protein networks elsewhere in the cell or cellular network. As cell-type specific interaction networks come on line in coming years, we look forward to integrating with our analyses.
Next, we tested the hypothesis that genes on 16p are over-represented in analysis of differential expression in the brains of individuals with ASD vs. controls. We identified differentially expressed genes between ASD cases (n = 51) and controls (n = 936) from a recent publication and retained those significantly variable at a bonferroni-significant level (n = 83 genes) (Gandal et al 2018 Science). We used a chi-squared test for over-representation of genes on 16p (n = 383 in Gandal dataset) in this n = 83 differentially expressed gene set. We did not find over-or under-representation of 16p-related genes in the Gandal DEG set (p > 0.05). Given the genetic heterogeneity of ASD, among the other non-genetic factors contributing to expression variability between ASD cases and controls, we do not find it surprising that there is no overlap between these gene sets.
In summary, we believe it is most likely that the genes on 16p -modulated by both common variants on 16p and the 16p11.2 deletion -are integrated in a complex network that is not ascertainable through canonical gene set enrichment approaches. We are engaging with members of the community to develop approaches to extract additional biological meaning out of regional variation in gene expression. We have added these analyses to the main text (page 6) and supplement (Supplementary Table 1 We agree this important section of the methods deserves additional detail in the text; we have expanded this methods section in the manuscript (page 18). In brief, for creating genomic blocks of 2,000 SNPs, we identified the first PRS SNP on chromosome 1 (the SNP closest to the first base pair), counted 2,000 PRS SNP, and called that the first partition. Then, we counted the next 2,000 PRS SNPs on chromosome 1, called that the next partition, etc, until we ran out of SNPs on chromosome 1. Then we started the same process on chromosome 2, etc. We repeated this for blocks of different sizes (3,000 SNPs, 4,000 SNPs, 5,000 SNPs, and 6,000 SNPs), as well as repeated the entire process starting at the ends of chromosomes and going backwards. The partitions were not based on LD blocks, nor on genes/gene density. Comment #3.3: " Figure 1E. The number of blocks removed stops at n=25. Is that just because there was no more effect?
In Figure 1E, the number of blocks removed stops at n = 25 because that is the number of these blocks located within 16p (median length of each block: 1.31 Mb). We have clarified this point in the text (page 19). Comment #3.4: "One would expect that the SE of the pTDT would get larger as more blocks are removed, but that doesn't appear to be the case." The SE of the S-pTDT decreases with larger sample size (right plot) but does not vary with the number of SNPs in the PRS partition (left plot). Comment #3.5: ( Figure 1E) "It also seems like there is a trend that may become significantly protective at one point. What would happen beyond =25? In other words, once authors remove the most over transmitted blocks, are there protective blocks in the 16p11.2 short arm?" Some of the blocks on 16p are (non-significantly) under-transmitted to ASD probands (Supplementary Figure  8), which indeed reflects a trend towards the common variants in that block being protective. This is not unique to 16p, but reflects a genome-wide pattern, where many regions of ASD common variation are under-transmitted in our three trio cohorts (Supplementary Figure 4). This reflects the genetic variability among our trio cohorts, where it is only with some probability at a given locus that an ASD proband inherited the liability-increasing haplotype. Comment #3.6: "As a sensitivity analysis, the authors performed an analogous analysis using a cohort of ADHD trios and an external ADHD GWAS and they did not replicate the finding in ADHD. One could argue that ADHD may be the worst condition to perform such a sensitivity analysis since the PRS doesn't explain much variance. Schizophrenia would appear to be much more relevant. The PRS is more robust, cohorts are larger, and the 16p11.2 locus is associated with schizophrenia." We agree that the SCZ PRS is a more predictive instrument for SCZ than is the ADHD PRS for ADHD. However, relative to the ASD PRS, the ADHD PRS performs well (ADHD Nagelkerke's R 2 = 5.5%, Demontis et al. 2019 Nature Genetics, vs. ASD Nagelkerke's R 2 = 2.5%, Grove et al. 2019 Nature Genetics). Regarding cohorts, the schizophrenia cohorts in the PGC are case-control design and not trio, which is required for the within-family transmission analysis of S-pTDT. In contrast, we were able to use ADHD trios from the PGC. Finally, the 16p11.2 locus is also associated with ADHD (Niarchou et al. 2019 Translational Psychiatry).

Comment #3.7: "The relationship between gene density and over-transmission is an important point and should be represented in a figure in the main text."
We agree this is an important point and have moved one of the panels relating gene-density and over-transmission to the main text as an inset to Figure 1F. Comment #3.8: "The authors asked whether, on average, the 200 neuronally expressed genes on 16p were differentially expressed in response to the 16p11.2 deletion. Genes on 16p had significantly lower expression in the deletion lines. The deletion's effect on 16p genes differed from the effect on all other 8,533 neuronally expressed genes in the genome (P = 0.02 -somewhat of a trend -), whose expression was not, on average, changed by the deletion (P = 0.43). Can the authors provide information on the non-neuronally expressed genes on the short arm? Would these represent a better "control group", providing a stronger contrast? Shouldn't the Y-axis of figure 2C represent "fold change" instead of t-stat?" We thank the reviewer for this thoughtful comment, and agree that the low expression condition should be included in the analysis for comparison. We have assessed the effect in low expression genes in both the isogenic 16p11.2 deletion lines, and using the regional PRS, and confirmed that the effect on decreased expression is attenuated in those lower expressed genes in both sets of analyses. We have summarized the findings in a figure below in this response. In the main text and supplement, we have added text describing the analyses and results for both the isogenic deletion (main text page 9, supplementary figure 16) and for the regional PRS approach (main text page 11, supplementary figure 20). Figure 2C, the left panel is in log(fold-change) for intuitive interpretability, while the right panel is the statistical comparison of the two groups that incorporates uncertainty in the fold-change estimates, hence displaying the changes in uncertainty-normalized t-statistics. Comment #3.9: "Increased ASD PRS within 16p was associated with decreased expression in glutamatergic neurons of genes through the 16p region. Do the authors observe an increase in the variance of gene expression? In other words, is this a simple shift in mean expression, or do results also suggest that there may also be some genes with an increase in expression?"

Regarding units in
Thank you for this interesting question. To explore this further, for each 33Mb region of the genome ("partition"), we associated regional PRS with expression of each gene and extracted the association t-statistic across our 544 samples. There was no association between either the partition's mean(t-statistic) or |mean(t-statistic)| and variance(t-statistic) across partitions (see plots below, where each dot is a partition with 16p in blue. P > 0.05 for each). For 16p specifically, we do not see a dramatic increase in expression variance given the decrease in expression averaged across all genes. Comment #3.10: "Authors show that the CNV-telomeric contacts (n = 291 100kb x 100kb contacts) are 2.9x more frequent than contacts between distance-matched control regions on 16p (n = 1,808 100kb x 100kb contacts, P < 1e-10). However, the 16p11.2 region is 30MB away from the telomeric region. I don't see, therefore how they can test distance-matched regions on the short arm. Distance-matched regions would only be found downstream of the 16p11.2 locus on the long arm." We define the telomeric region from 0 Mb to 5.2 Mb on chromosome 16, while the 16p11.2 (proximal) CNV ranges from 29.5Mb-30.2Mb. The contacts between these regions are denoted in Figure 4C in the blue shaded rectangle inside the larger triangular contact matrix. The range of distances encompassed between these contacts begins at 24.3Mb in distance (contact between Mb 5.2 of the telomeric region and Mb 29.5 of the CNV: 29.5 -5.2 = 24.3) and extends to 30.2Mb in distance (contact between 30.2Mb of the CNV to 0 Mb of the telomeric region). Therefore, the distance-matched control regions are contacts on 16p that span 24.3Mb to 30.2Mb. There are many such 100kb x 100kb contacts (n = 1,808), and such contacts are denoted in the red shaded trapezoid in Figure 4C (for example, contact between Mb 6 and Mb 31 = 25 Mb apart). Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Decision Letter, first revision:
Our ref: NG-A59672R 4th August 2022 Dear Dan, Your revised manuscript "Statistical and functional convergence of common and rare genetic influences on autism at chromosome 16p" (NG-A59672R) has been seen by the original referees. As you will see from their comments below, they find that the paper has improved in revision, and therefore we will be happy in principle to publish it in Nature Genetics as an Article pending final revisions to address the referees' remaining points and to comply with our editorial and formatting guidelines.
We are now performing detailed checks on your paper and we will send you a checklist detailing our editorial and formatting requirements soon. Please do not upload the final materials or make any revisions until you receive this additional information from us.
Thank you again for your interest in Nature Genetics. Please do not hesitate to contact me if you have any questions.
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reviewer #1 (Remarks to the Author):
Reviewer #1: Comment #1.1: "The deletion described and studied in this manuscript is the 16p11.2 proximal deletion explaining up to 1% of autistic cases. It is particularly interesting that distal to this locus is another recurrent deletion, 16p11.2 distal deletion, not mentioned in the manuscript. That deletion, the 16p11.2 distal deletion confers high-risk of very similar phenotypes (including autism, cognitive impairments and obesity). If that deletion, the 16p11.2 distal deletion, also affects gene expression on 16p in similar manner as the 16p11.2 proximal deletion and the 16p ASD PRS, the story would be more convincing. Thus, my recommendation is to include as well data on the 16p11.2 distal deletion in the manuscript."

Authors reply
We thank the reviewer for this thoughtful comment about the possibility that other ASD-associated CNVs -especially those located on 16p -may confer disease liability through similar mechanisms to the proximal 16p11.2 studied here. The reviewer notes the potential relevance of the 16p11.2 distal deletion, and we agree this proximate and disease-associated CNV would be interesting to investigate. However, as far as we can tell, there are no published whole-genome RNA-sequencing datasets of the 16p11.2 distal deletion, including in the provided reference of Sønderby et al. It is also materially infeasible for us to generate additional isogenic distal deletions at this time. Without RNA-sequencing, we are unable to evaluate our model for this deletion. That said, we are very interested in whether other neuropsychiatric CNVs confer disease liability through similar mechanisms to those described for the proximal 16p11.2 CNV. We are actively testing this hypothesis in our research group, and we look forward to sharing our findings with the community in the future. Finally, we clarified in the manuscript that we analyzed the proximal and not the distal deletion (page 3).
Further comments from Reviewer #1: I find the results presented in this manuscript most interesting. However, they should be confirmed. Authors can identify RNA-sequenced samples suitable for confirming their findings or they can analyze the isogenic neuronal cell lines with CRISPR/Cas9-mediated "distal" deletion of 16p11.2. That may reveal that the deletion also associates with depressed average gene expression across 16p and, hence, confirm the findings.
Reviewer #3 (Remarks to the Author): The authors responded to all comments and questions in a satisfactory way.
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
The only response that remains unclear relates to comment 3.10. The authors describe the contacts between the 16p11.2 region and the telomeric region ranging from 24.3 to 30.2 Mb in distance. They give an example of the distance-matched control regions between MB6 and MB31. However, this contact is beyond the 16p11.2 region. Does this mean that the distance-matched control contacts do not include the 16p11.2 region?

Author Rebuttal to Initial comments
Dear Kyle, Thank you for sharing the reviewer comments. Please see our responses below:

Response to reviewer #1
We are glad the reviewer finds our manuscript of great interest. We interpret the reviewer's specific request here as asking for analysis of an additional isogenic CRISPR-generated lines of either the proximal or distal 16p11.2 deletion. While we've reviewed the literature and inquired broadly, such a resource does not seem to currently exist, unfortunately.
We're happy to reflect more on the one result in question (the data in Figure 2). We do find support for the observation from other analyses presented in the manuscript, including a) convergence with expression effects to the 16p ASD PRS ( Figure 4A), b) elevated chromatin contact between the 16p11.2 deletion region and the telomeric region of convergent effect ( Figure 4C), and c) lack of a similar observation at the 15q locus supporting the specificity of the 16p11.2 deletion effect on regional gene expression (Supplementary Figure 17). We will also add a note to the discussion that replication using further isogenic lines or very large patient-derived samples will be valuable once those resources are developed.
Should the NG editors or the reviewer have other questions or suggestions, we're happy to discuss.

Response to reviewer #3
Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
We are glad the reviewer finds our responses satisfactory. With regards to comment 3.10: yes, that is correct that almost all of the control contacts do not include the 0.7Mb 16p11.2 region. Three-dimensional contact frequencies as assayed by Hi-C are strongly dependent on the distance between the contact loci (decaying with distance). Thus, we defined control regions based on their distance in such a way that the range of control contact distances (red trapezoid in Figure 4C) is the same as the range of contact distances between the 16p11.2 region and the telomeric contact region (blue rectangle in Figure 4C). Hopefully this is clarifying.

Final Decision Letter:
In reply please quote: NG-A59672R1 Weiner 15th September 2022 Dear Dan, I am delighted to say that your manuscript "Statistical and functional convergence of common and rare genetic influences on autism at chromosome 16p" has been accepted for publication in an upcoming issue of Nature Genetics.
Over the next few weeks, your paper will be copyedited to ensure that it conforms to Nature Genetics style. Once your paper is typeset, you will receive an email with a link to choose the appropriate publishing options for your paper and our Author Services team will be in touch regarding any additional information that may be required.
After the grant of rights is completed, you will receive a link to your electronic proof via email with a request to make any corrections within 48 hours. If, when you receive your proof, you cannot meet this deadline, please inform us at rjsproduction@springernature.com immediately.
You will not receive your proofs until the publishing agreement has been received through our system.
Due to the importance of these deadlines, we ask that you please let us know now whether you will be difficult to contact over the next month. If this is the case, we ask you provide us with the contact Open Access This file is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. In the cases where the authors are anonymous, such as is the case for the reports of anonymous peer reviewers, author attribution should be to 'Anonymous Referee' followed by a clear attribution to the source work. The images or other third party material in this file are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. information (email, phone and fax) of someone who will be able to check the proofs on your behalf, and who will be available to address any last-minute problems.
Your paper will be published online after we receive your corrections and will appear in print in the next available issue. You can find out your date of online publication by contacting the Nature Press Office (press@nature.com) after sending your e-proof corrections. Now is the time to inform your Public Relations or Press Office about your paper, as they might be interested in promoting its publication. This will allow them time to prepare an accurate and satisfactory press release. Include your manuscript tracking number (NG-A59672R1) and the name of the journal, which they will need when they contact our Press Office.
Before your paper is published online, we will be distributing a press release to news organizations worldwide, which may very well include details of your work. We are happy for your institution or funding agency to prepare its own press release, but it must mention the embargo date and Nature Genetics. Our Press Office may contact you closer to the time of publication, but if you or your Press Office have any enquiries in the meantime, please contact press@nature.com.
Acceptance is conditional on the data in the manuscript not being published elsewhere, or announced in the print or electronic media, until the embargo/publication date. These restrictions are not intended to deter you from presenting your data at academic meetings and conferences, but any enquiries from the media about papers not yet scheduled for publication should be referred to us.
Please note that Nature Genetics is a Transformative Journal (TJ). Authors may publish their research with us through the traditional subscription access route or make their paper immediately open access through payment of an article-processing charge (APC). Authors will not be required to make a final decision about access to their article until it has been accepted.  Find out more about Transformative Journals Authors may need to take specific actions to achieve  compliance with funder and institutional open access mandates. If your research is supported by a funder that requires immediate open access (e.g. according to Plan S principles), then you should select the gold OA route, and we will direct you to the compliant route where possible. For authors selecting the subscription publication route, the journal's standard licensing terms will need


Identification of potential target genes in Homo sapiens, by miRNA of Triticum aestivum: A cross kingdom computational approach

Plant-derived miRNAs can be found in the human body after dietary intake, and they can affect post-transcriptional gene regulation in human. It is important to identify targets to determine the possible effects in human genes by using computational approach. In this study, 787 possible mRNAs human targets were predicted by 84 miRNAs of wheat. A total of 14 miRNAs were identified with individual binding to 33 mRNAs associated with schizophrenia, epilepsy, neurodevelopmental disorders, and various cancers, located in the 3′UTR of the mRNA. A functional enrichment was carried out, where the results showed associations to pathways such as dopaminergic synapse (hsa04728), and signaling pathways, significantly associated with the target genes. The prediction of target mRNAs in humans by wheat miRNAs, offer candidates that could facilitate the search and verification, which could be of relevance for future projects and therefor contribute in the therapeutic treatment of various human diseases.

Introduction
The microRNAs (miRNAs) of plants and animals are an abundant class of small non-coding RNAs sequences around ~22 nucleotides long [1]. They are derived from hairpin precursors, which have remained evolutionarily similar in various species, and have the ability to negatively regulate gene expression at the post-transcriptional level by complementary base binding in various regions of the mRNA (3 ′ UTR, 5 ′ UTR and ORF) and via the RNA-induced silencing complex (RISC) whose binding causes mRNA cleavage or translation inhibition [ [2][3][4]].
It is expected that 30-60% of the genes in humans are miRNA regulatory targets [5], involved in various developmental processes, molecular, cellular pathways and human disease [6]. When crossing the intestinal barrier after ingestion miRNAs can be packaged in microvesicles, circulate through the human body and be identified in fluids such as human plasma [7] and breast milk [8].
Identifying mammalian gene targets in a timely manner by plant miRNA offers a reliable, effective and economical way to determine the possible effects of interactions given by exogenous molecules, regardless of their origin, on human mRNA [ [9][10][11][12][13][14]].
Recently, plant miRNAs and their interactions have been actively studied, contributing to the knowledge of coevolution mechanisms between plant miRNAs and human mRNAs. For example, the in silico prediction on the possible regulation of plant miRNA on human mRNA was tested, coincided with an increase in cholesterol-LDL levels in vivo and in vitro, caused by the decrease of LDLRAP1 protein expression, providing the first evidence that food miRNA could bind to a target mRNA and influence the regulation of interspecies expression [15]. Similarly, through bioinformatics analysis, found a number of potential target genes for MIR2911 from the Lonicera japonica plant (honeysuckle) were found in different viruses such as H1N1, H3N2, H5N1 and H7N9. Later, in cultures of cellular lines infected with H1N1 and transfected with synthetic MIR2911, they managed to significantly decrease the viral count because the target genes: PB2 and NS1, are essential for replication. Subsequent functional studies showed that miR2911 inhibited several IAVs (influenza A viruses) including H1N1, H5N1 and H7N9, decreasing mortality in an animal model [16]. Suggesting the first reports that exogenous miRNAs regulate gene expression in endogenous mRNAs similar to mammalian miRNAs and in turn the use of the bioinformatics approach, playing an important role in target identification.
In the present study we have selected the monocotyledonous plant Triticum aestivum (wheat) because it is a common source in human nutrition, we performed an in silico prediction analysis of potential mRNAs target in Homo sapiens (human) by mature miRNAs of wheat, and its relationship to diseases through an enrichment of the possible target genes. In addition, we identified target genes with diverse functions in cellular processes, providing information on regulatory mechanisms given by wheat miRNAs. Also, a functional enrichment has been carried out to explore the participation of the target human mRNAs involved in biological processes and associated pathways, including apelin signaling pathway, fanconi anemia pathway and dopaminergic synapse.

Data collection
The complete set of known mature miRNA sequences, which includes a total of 122 miRNAs from wheat, were downloaded from the registered database (version 22.1) (http://www.mirbase.org) in FASTA format. To reduce errors in prediction and future detection, we applied rigorous filtering for mature miRNAs, i) The mature miRNA should have more than 1400 reads during sequencing reported in its bioinformatic analysis; ii) The miRNA size should be in the range of 19-25 nucleotides and iii) All mature miRNAs should be in the same direction from 5 ′ to 3'. Following this, 84 miRNAs from wheat were selected as reference data set for the study.

Prediction of potential targets of miRNAs
The prediction of possible targets on human was carried out by the sequence similarity search tool BLASTn (2.9.0) on its NCBI website (http://www.ncbi.nlm.nih.gov) using the default settings, an E-value threshold of 10; and the search coordinates were set in the nucleotide sequence 2 to 13 of each miRNA. From the RefSeq RNA database the predicted miRNA target sites were identified by aligning the sequences in a complementary way to the miRNA of wheat with human mRNA. RNAhybrid [17] was used to scan the hybridization patterns of the individual miRNAs/mRNA following Zhou Z. et al. methodology with some modifications [16]. i) The seed region had to perfectly match the nucleotide 2 to 7 of the mature miRNA. ii) The minimum free energy (MFE) for the hybridization was below − 25 kcal/mol iii) No G-U bonds in the seed region. We also used the MIR168a as a target prediction control [15].

miRNA homology
We searched for homology using MirCompare between 2042 H. sapiens and 84 T. aestivum matured miRNAs from miRBase was used [18]. The parameters used for this analysis were taken from Minutolo et al., 2018 [19], with modification in the threshold of the seed region according to the guidelines provided by the default settings of the tool. The possible target genes of the resulting human miRNAs with homology to wheat miRNAs as reported by MirCompare were found by an analysis with TargetScan 8.0 algorithms (http://www.targetscan.org).

Functional enrichment analysis
The functional analysis of the targets, including biological processes was done using clusterProfiler [20], package inside programming language R version 4.0.3 [21]; with only those target sequences with miRNA binding sites in the 5 ′ UTR, 3 ′ UTR and CDS regions. Designating the potential processes with an adjusted p-value lower than 0.05 (FDR <0.05) that involve the predicted mRNA sequences and suggesting the probable influence in the regulation of plant miRNAs in the human organism.

Characteristics of the binding of wheat mature miRNAs with human mRNAs
According to the methodology previously mentioned, the results obtained were similar to the ones described by Zhang et al. 2012 [15]. We selected 84 mature miRNAs from wheat and predict putative targets genes in humans, based on the criteria specified in methods section, which resulted in a total 889 miRNA/target interactions with 787 genes with diverse functions, the characteristics and specific target sequences of the 84 miRNAs analyzed in this study are provided in the (Supplementary Table 1). 70 of these genes had interactions with more than one miRNA and 6 come from miRNAs families, as previously reported by Han R. et al. [22] (Table 1). These miRNAs were divided in two groups; in group one, 14 miRNAs identified with one to four target genes, giving a total of 30 target genes. In group two 70 miRNAs with five or more target genes, giving a total of 757 target genes. The tae-miR9652-3p and tae-miR9652-5p of group one, originated from the same precursor (pre-miRNA) in miRNA biogenesis, have binding sites in the mRNAs of different genes. In the second group the tae-miR9657b-3p and tae-miR9657b5-p, were found as a pair originating from the same pre-miRNA, obtaining 18 targets for miR-3p, and 10 for miR-5p strand, all being different. The miRNAs with the highest number of targets were tae-miR9657c (19 genes), tae-miR9654a (21 genes), tae-miR160 (22 genes), tae-miR396 (23 genes), tae-miR9662a (23 genes), tae-miR9666a (24 genes), tae-miR9664 (28 genes), tae-miR9655 (31 genes), tae-miR9666b-3p (34 genes), tae-miR164 (35 genes).
In order to show the targets and MFE for the binding between some wheat miRNAs with human mRNAs, data was filtered to show the targets genes with lower value of MFE and were the sites union are located in 3 ′ UTR binding ( Table 2). 10 of the miRNAs show 1 target gene, 2 miRNAs show two targets, and 2 miRNAs show three and six targets, with a MFE between − 34.1 and − 43 kcal/mol. The tae-miR160 shows 10 target genes with a free energy value between − 36.9 and − 45.9 kcal/ mol.

Analysis of human homologs of wheat miRNAs
In the analysis we found homology candidates in 9 human miRNA sequences, with 7 wheat miRNA sequences. Supporting this, we identified target gene predictions of human miRNAs using the TargetScan tool and determined those targets in common predicted in this study as a whole (Table 3).

Functional enrichment of target genes by wheat miRNAs
Functional annotation of target genes is an important step in understanding the respective functions. A gene ontology (GO) enrichment analysis was performed in the Kyoto Encyclopedia of Genes and Genomes (KEGG) to identify the functions that miRNAs could have on their predicted targets in humans, where all mRNAs that had binding sites in the prediction were used. From the complete target search data, 670 different genes were used, selected according to probable binding site, including only those in the 5 ′ UTR, 3 ′ UTR, and CDS region and excluding those in non-coding sequences for functional enrichment analysis. The main GO terms classified in biological process are shown, with significance adjusted p-values (<0.05) indicating that a total of 17 mRNAs enriched the biological process of homophilic cell adhesion (GO:0007156), through adhesion molecules to the plasma membrane, and 22 mRNAs enriched the biological process of cell-cell adhesion through adhesion molecules to the plasma membrane (GO:0098742), by miRNAs known from wheat in the prediction (Fig. 2).

Discussion
Based on the results of this work and the extensive study on the mechanisms of plant miRNAs when entering the circulatory system through the diet, entering the cells to achieve a target and regulate the expression of a gene [15]. Some authors [23,24], suggest the use of available and evaluated tools for predictive analyses, since they recognize the fundamental importance of identifying the target genes in the human genome by exogenous miRNAs, in order to infer their functions and regulatory mechanisms, however they agree that it is convenient to explore and combine predictive tools because individually, they may not reflect the abundance and possible regulatory capacity of miRNAs, improving the speed and accuracy of the targets. In our research we used two search engines, that way we were able to simulate a possible interaction between species. By using BLASTn as the first instance as a search engine for complementary sequences with the updated version of the RefSeq human transcriptome database, we were able to find several targets and their variants, which is important because factors related to expression variation in target genes can provide a solid foundation and more systematic view of the genetic mechanisms that are related to human genetic variation [25].
The results obtained by the BLASTn analysis were filtered in a second computational method using RNAhybrid tool as used successfully used in several studies of cross kingdom miRNAs [16,26,27], which allows us to determine the MFE, characteristics and complementary regions of the miRNA/mRNA interaction, which allowed filtering of those interactions that we consider energetically favorable (MFE-25 kcal/mol). This provided us with a combination of algorithms in the study that allowed us to obtain a total of 787 genes and their variants with different functions as shown in (Supplementary Table 1). In addition, Table S1 shows the transcripts of non-coding RNAs (ncRNAs) that were targeted by the 84 miRNAs, since miRNA targets are not limited to coding transcripts, but can also be found in ncRNAs, indicating that miRNA targeting mechanisms can be complex by binding to a wide variety of targets [28].
A total of target genes for 84 miRNA of wheat, divided into two groups according to the number of targets, and families that each miRNA obtained, since miRNAs with a larger number of target genes should be monitored in human fluids, due to their possible ability to significantly regulate human receptor cells [9]. These targets were located in various regions along the mRNA (5 ′ UTR, CDS, 3 ′ UTR), increasing the area of study on the interactions given by exogenous molecules, it could help understand the role of binding sites in the mRNA with the miRNAs [29]. We also evaluated the type of binding, obtaining 100% of the possible binding as imperfect base-pair by at least 1 base, however the analysis was made with the criterion of perfect binding of the seed region which would suggest that the possible mechanism of regulation of expression according to the reported silencing by mature miRNAs [2] could be translation inhibition and not degrading by cleavage the mRNA, which is sufficient for the miRNAs to suppress their targets [23,30]. In plants, miRNA bind to their targets with high complementarity, causing messenger degradation, however in animals, miRNA recognition in just their seed region is extremely important in target recognition, this complementarity can trigger expression regulation, [31,32]. Even in humans, the seed region of some miRNA can define toxicity or tumor-suppressive for some cancer cell lines [33]. This is due to seed region target recognition mechanisms [34,35], show that an initial search for target sites transitions from being transient to kinetically stable, when target complementarity binds to the full length of the seed region of miRNA [36]. Furthermore, several authors have identified and demonstrated regulation of target gene expression, even when miRNAs have not bound in a fully complementary manner, showing evidence of their function [16,37,38].
In this study we included a homology analysis, identifying 7 wheatderived miRNAs show homology with 9 human miRNAs, revealing a series of target genes in common, through our prediction. Although plant sequences do not show high homology with human sequences, results obtained by author [39], show that they can perform similar functional activities to endogenous sequences, in a cross-kingdom interaction [19]. Therefore, we consider this dataset provides a solid basis for future projects.
In our mRNA and miRNA interaction analysis we aimed to identify the genes that are the target of multiple miRNAs and their possible health effect through KEGG pathways. We found that the target genes are associated with the calcium signaling pathway, showing participation of the CD38 gene, associated with diseases such as asthma, neoplasia and neuroimmune diseases [40], as well as the ADCY1 gene, involved in insulin secretion, related to diabetes mellitus and its correlation with pancreatic cancer [41].
Among the pathways, three stand out, the pathway of fanconi anemia (hsa03460), circadian entrainment (hsa04713), and dopaminergic synapse (hsa04728), which we have found a relation with a couple of publications [42,43](Gluten-Free Diet for the Treatment of ADHD; Pilot Study) where in a case study they describe the difficulty of finding a diagnosis to a series of neuropsychiatric manifestations, described as "non-celiac gluten sensitivity" (NCGS), where a series of symptoms related to the consumption of gluten-containing foods (wheat) are characterized in people not affected by celiac disease (CD) or wheat allergy (AW).). Attention deficit hyperactivity disorder (ADHD) related to the circadian entrainment pathway [44], and schizophrenia related to the dopaminergic synaptic pathway, are some of the main symptoms related to wheat consumption in case studies, in addition, including a low-ferritin anemia in the clinical picture. Despite attempts to consolidate a possible explanation, this has not been achieved due to lack of information on response mechanisms or genetic components of the disease and therefore the number of people who may be affected is uncertain. In the enrichment results, we have 5 genes, GNB5, GNAI2, PRKACB, CALM1, and CAMK2D as possible targets in common for miRNA regulation of wheat involved in the dopaminergic synapse pathway and circadian entrainment.
Since the report of miRNA plants, regulating gene expression by miR168a in mice and humans [15], the identification of new targets by miRNA is of major importance, due to its high efficiency, low cost, and time saving [16,38,45]. Since this could be a way to explain negative events in human health, related to some diseases that could risk public health, such as cancer, schizophrenia, and metabolic diseases.

Conclusion
Considering that wheat and humans are different species, it is difficult to determine the possible alterations given by exogenous molecules caused by cross-kingdom interactions, however, it is necessary to highlight the importance of the use of bioinformatics tools, for the prediction of possible targets, because this study demonstrates the possible effects of wheat miRNAs in regulating complex disease networks like a schizophrenia, epilepsy, tumor suppressor, inflammation, and diverse type of cancer. Due to the fact that the information collected could facilitate the search and verification, that could be of importance for future projects. In addition, these findings can generate new lines of research in various areas, such as nutrition, medicine, biotechnology, agriculture, etc. It can also show possible miRNAs candidates that can further contribute in therapeutic treatment of various human diseases.

Declaration of competing interest
None.