text:
Imputation of Unordered Markers and the Impact on Genomic Selection Accuracy

Genomic selection, a breeding method that promises to accelerate rates of genetic gain, requires dense, genome-wide marker data. Genotyping-by-sequencing can generate a large number of de novo markers. However, without a reference genome, these markers are unordered and typically have a large proportion of missing data. Because marker imputation algorithms were developed for species with a reference genome, algorithms suited for unordered markers have not been rigorously evaluated. Using four empirical datasets, we evaluate and characterize four such imputation methods, referred to as k-nearest neighbors, singular value decomposition, random forest regression, and expectation maximization imputation, in terms of their imputation accuracies and the factors affecting accuracy. The effect of imputation method on the genomic selection accuracy is assessed in comparison with mean imputation. The effect of excluding markers with a large proportion of missing data on the genomic selection accuracy is also examined. Our results show that imputation of unordered markers can be accurate, especially when linkage disequilibrium between markers is high and genotyped individuals are related. Of the methods evaluated, random forest regression imputation produced superior accuracy. In comparison with mean imputation, all four imputation methods we evaluated led to greater genomic selection accuracies when the level of missing data was high. Including rather than excluding markers with a large proportion of missing data nearly always led to greater GS accuracies. We conclude that high levels of missing data in dense marker sets is not a major obstacle for genomic selection, even when marker order is not known.

some loci to have zero coverage in some individuals (Elshire et al. 2011). The proportion of missing data depends on the sequencing depth and library complexity. Greater sequencing depth leads to a smaller proportion of missing data but increases genotyping cost. Less-complex libraries, on the other hand, will have less missing data but a fewer markers. To generate a large number of markers at low cost, low sequencing depth is commonly used, leading to a large proportion of missing data points. Most analyses require a complete dataset; therefore, marker imputation is a necessary step before GBS data can be used for most purposes.
Imputation has been shown to increase power in association mapping studies (Marchini et al. 2007;Marchini and Howie 2010) and, for GS, imputation can enable the use of low-density genotyping without a major loss in accuracy because a closely related reference panel genotyped at high density can be used to impute markers not present in the low-density marker panel. (Habier et al. 2009;Weigel et al. 2010;Dassonneville et al. 2011;Mulder et al. 2012). Although several highly accurate and widely used imputation algorithms have been developed to assign allelic states of missing values in genotype data (reviewed by Pei et al. 2008 andMarchini et al. 2010), these algorithms were designed for human genetic data, and they require that the order of the markers be known because they are based on constructing haplotypes. For species lacking a reference genome and complete reference linkage map such as wheat, Triticum aestivum L., the majority of markers typed on a given population are unordered, and current genotype imputation methods cannot be used. Although for biparental populations linkage maps can be constructed, breeding populations for GS are derived from multiple parents and not well structured for developing genetic maps. Thus, alternative imputation strategies that are map-independent are necessary when GBS is used for species lacking a reference genome sequence and for populations unsuitable for linkage map construction. There are many general imputation methods that do not require any prior information about the variables to be imputed. Although these methods are used across many disciplines, they have not been tested for imputation accuracy of genome-wide marker data. It is also not known how imputation with a general and potentially less accurate method before GS model training will affect the GS model accuracy. However, we expect these imputation methods to improve the GS accuracy because during the imputation step, genotypic information from both the training and selection sets is used to estimate missing values. Thus, the validation set helps improve imputation of the training set and vice versa.
The objective of this study was to evaluate imputation strategies that do not require previous information about the order of the markers. The imputation methods compared were: mean imputation (MNI), k-nearest neighbors imputation (kNNI) (Troyanskaya et al. 2001), singular value decomposition imputation (SVDI) (Troyanskaya et al. 2001), expectation maximization imputation (EMI) (Dempster et al. 1977), and random forest regression imputation (RFI) (Stekhoven and Bühlmann 2011). By using array-based genotypic datasets with varying levels of simulated missing data, we compared these methods in terms of their imputation accuracy, computation time, and impact on GS prediction accuracy. The factors affecting imputation accuracy for each method at the marker genotype and individual genotype level were also examined. Finally, we determined whether excluding rather than including markers with high levels of missing data could lead to higher accuracy.

Original datasets
We used five different datasets consisting of genome-wide markers and breeding value estimates. These datasets are referred to as winter wheat (WW), spring wheat (SW), drought tolerant maize (DTM), North American barley (NAB), and stem rust resistant wheat (SRRW). The WW data (Supporting Information, File S2) consists of 374 elite inbred individuals originating from the Cornell winter wheat breeding program. The markers consisted of 1158 polymorphic diversity array technology (DArT) (Akbari et al. 2006) markers coded as "21" and "1." For a more detailed description of this dataset, refer to Heffner et al. (2011). The traits used for the evaluation of crossvalidated GS accuracies for WW were grain yield, height (HT), protein, and days to heading. The SW data are a historical dataset consisting of 599 elite inbred spring wheat lines originating from the International Maize and Wheat Improvement Center (CIMMYT) wheat breeding program. The markers consist of 1279 polymorphic DArT markers coded as "0" and "1" and the trait used for the evaluation of cross-validated GS accuracies was grain yield in CIMMYT megaenvironment 1. The DTM data consists of 264 tropical CIMMYT maize lines. The trait used to calculate cross-validated GS model accuracies for DTM was grain yield. The marker data consists of 1135 single-nucleotide polymorphisms coded as "21", "0", and "1." For more details about the SW and DTM datasets, or to access these datasets, refer to Crossa et al. (2010) The SRRW dataset consists of 360 recent, elite CIMMYT spring wheat lines that have been selected for quantitative resistance to stem rust caused by Puccinia graminis f.sp. tritici. The markers consist of more than 130,000 GBS polymorphisms. Three different versions of the SRRW GBS data, described in Table 1, were created based on different per-marker percent missing data thresholds. For the first version referred to as SRRW version NA20 (File S3), markers were excluded if they had more than 20% missing values, which resulted in 2014 total markers. For the second set and third sets, referred to as SRRW versions NA50 and NA70 (File S4 and File S5), markers were excluded if they had more than 50% and 70% missing data, respectively, and then 2014 markers were randomly selected. The percent of the data points that were missing in the original WW, SW, DTM, and NAB datasets was between 0.2 and 3%. This low level of pre-existing missing data was assumed to have a negligible effect on the imputation and GS accuracies and for these datasets the original marker data are referred to as version NA0.

Calculation of LD between marker pairs
For the original WW, SW, DTM, and NAB datasets, LD between all marker pairs was measured using the r 2 statistic, where r 2 between two markers was calculated using the formula: where D = x 11 -p 1 p 2 ; x 11 is the probability of observing the combination of allele 1 at marker j and allele 1 at marker l, p 1 is the probability of allele 1 at marker j, q 1 is the probability of allele 2 at marker j, p 2 is the probability of allele 1 at marker l, and q 2 is the probability of allele 2 at marker l. A maximum likelihood estimate of x 11 was obtained using an expectation maximum approach reviewed by Foulkes (2009). All calculations of the r 2 statistic were implemented in the R package "genetics" (Warnes et al. 2011).

Missing data simulation
For each of the WW, SW, DTM, and NAB datasets three versions of the genotypic data, summarized in Table 1, were created with different levels of simulated missing data. In each of the versions: NA20, NA50, and NA70, missing values were introduced at random but the maximum percent missing data at a given marker was set to 20%, 50%, and 70% respectively. Examples of the simulated markers sets are illustrated in the Figure S1. A total of 10 replicates of each simulated dataset were created, and the mean percent of total data points that are missing across the 10 replicates is shown in Table 1. The distribution of per-marker percent missing values from the SRRW data versions NA20, NA50, and NA70 were used to assign the percent missing at each marker for each of the WW, SW, DTM, and NAB datasets to produce versions NA20, NA50, and NA70, respectively. Across all the missing data versions of all the datasets, the percent missing per marker distribution had a long left tail and a large concentration of values near the threshold level.

Imputation methods
In all cases, the genotypic data were considered continuous variables. The methods MNI, kNNI, SVDI, EMI, and RFI were used to impute the simulated missing values. For all methods the input was an m x n genotype matrix M with m individuals and n markers. For MNI, each missing data-point x ij at a given marker j was replaced with the mean of the non-missing values at that marker.
For kNNI (Troyanskaya et al. 2001), the data points were imputed by replacing them with the weighted average of the data points at the k closest markers. Euclidean distance was used as the measure of marker distance. Euclidean distance between marker genotype vectors q ! and v ! of length m was defined as: In detail, (1) missing values were first replaced using MNI and the Euclidean distance between all of possible pairs of marker vectors was computed. Each marker was included in the marker matrix twice, both in its original and flipped state, to ensure that markers in negative LD would not be considered distant to the marker of interest.
(2) For each marker j, markers were sorted based on Euclidean distance to marker j.
(3) For each row i of marker j the weighted average of the k closest markers with nonmissing values at row i were used as an estimate of marker data point x ij . The weight of each marker was 1/d 2 were d is the Euclidean distance between marker j and the marker to be weighted. kNNI makes no assumptions about the distribution of the data.
For SVDI (Troyanskaya et al. 2001), a singular value decomposition of genotype matrix M was used to obtain a set of the k most significant Eigen-vectors of the markers. These k Eigen-vectors were then used as the predictors for linear regression estimation of the missing data points. SVDI was implemented in R (R Development Core Team 2011) using the package "bcv" (Perry 2009). The genotype matrix M can be described as: Where U has dimensions m · k, V has dimensions n · k, and S is a k · k diagonal matrix. U contains the left singular vectors with are equivalent to the Eigen-vectors of the markers. The corresponding singular values are in the diagonal elements of S. The singular values are equivalent to the square root of the Eigen-values. The k most significant Eigen-vectors of the markers were those with the k largest Eigen-values. The imputation procedure is described as follows: (1) Missing values were originally imputed using MNI.
(2) Singular value decomposition was used to estimate the k most significant Eigen-vectors of the markers:Û.
(3) For each marker j, linear regression coefficients of each column ofÛ were estimated by the multiple linear regression equation: Where Y is a column vector for marker j,Û is an m · k matrix of k Eigenvectors, b is a vector of regression coefficients and e is a random error term. Only individuals with nonmissing values in Y were used to estimate b. (4)Û and the estimates of the regression coefficients,b, were used to estimate the missing values at marker j. (5) Using the current version of the genotype matrix, we repeated steps two through four for a total of 10 iterations, sufficient to meet the convergence criteria, which was: RSS is the residual sum of squares between the nonmissing values and their SVDI model approximation. RSS 0 and RSS 1 are the RSS values of successive iterations. SVDI assumes that the genotype matrix is multivariate normal distributed. For the optimal k value calculation methods and results for both kNNI and SVDI, see File S1.
Optimal k values are listed in Table S1. For EMI, the nonmissing marker data were used to obtain maximum likelihood estimates of the vector of means,û, and covariance matrixX of the individuals based on the markers. These estimates were then used to obtain multiple linear regression estimates of the missing marker values.û, andX were then re-estimated and were used to re-estimate the missing marker values. This process was repeated until the difference between the new estimate and the previous estimate ofû þX ÃX T was 0.02 or less. EMI was implemented using the R package rrBLUP (Endelman 2011). For a more detailed description of this EMI algorithm, refer to . For a more through description of the EM imputation algorithm in general, refer to Dempster et al. (1977). For RFI, missing marker values were estimated using random forest regression (Breiman 2001) using all available data to predict the missing values for every marker. RFI was implemented in R using the package "MissForest" (Stekhoven and Bühlmann 2011). The imputation procedure was: (1) for marker matrix M, markers were sorted from lowest to highest percent missing and missing values were imputed using MNI. (2) At each marker j containing missing values, the nonmissing values, Y, were used to grow 100 random forest regression trees Q1. . . Q100. Each tree was grown using a bootstrapped sample of individuals Y and a random sample of ffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 1 p marker predictors were used where n-1 is the number of markers excluding marker j. Each tree Q contains the terminal node values and a set of instructions for recursively partitioning the observations into the terminal nodes: these instructions include the split variables at each node, and the value of the split variable used for partitioning. (3) Missing values at marker j were imputed as: where x is an input vector.
(4) Marker j was then updated in marker matrix M by using theŶ values as the estimate of the missing values. (5) Steps two through four were repeated for each subsequent marker until all markers were imputed. (6) Then, using this imputed matrix, we repeated steps two through five until convergence or for a maximum of 10 iterations. Convergence was declared as soon as the DN increased for the first time where: M 1 and M 0 are the newly imputed and previously imputed marker matrices respectively. If the convergence criterion was met, M 0 was used as the final estimate of M. RFI makes no assumptions about the distribution of the data. The implementation of all imputation methods is demonstrated in File S6.

Imputation accuracy calculations
The per-marker imputation accuracy, R 2 m , was described using the R 2 value between predicted data points and the original data points for a given marker vector or individual vector x of length j. The R 2 was defined as The R 2 m , as well as the imputation R 2 of the individual genotypes, referred to as R 2 i , were calculated. For each dataset and missing data level, average R 2 i and R 2 m across the 10 missing data simulations were also calculated and referred to as R 2 i and R 2 m . To compare with imputation accuracies reported in other publications, for each R 2 m value, the equivalent percent correct was also calculated. Because imputed values were continuous, the percent correct for each marker could not be directly calculated. Instead, for each marker, equivalent percent correct values were determined by simulation using each marker's MAF and R 2 m (see File S1).

Computational time
For the first replicate of simulated missing datasets, whenever a dataset was imputed, the number of seconds required for imputation to be completed using one central processing unit was recorded. All jobs were submitted to the Computational Biology Service Unit at Cornell University, which uses (1)

Assessment of factors affecting imputation accuracy
For each imputation method factors affecting the imputation accuracy were assessed. A marker's minor allele frequency (MAF), number of nonmissing data points, and level of LD with other markers were considered as factors that could impact its imputation accuracy. The distance between an individual and its closest relative and the expected prediction error variance (PEV) were considered as factors affecting the imputation accuracy on an individual genotype basis. The impact of each of these factors was assessed for each imputation method using the WW, SW, DTM, and NAB datasets post imputation. First, the impact of MAF on the imputation accuracy was assessed. For each dataset-imputation method combination, R 2 m was averaged across dataset versions NA20, NA50, and NA70 and this overall estimate of marker imputation accuracy is referred to as R 2 m .The median R 2 m for each value of MAF rounded to the nearest tenth was calculated. The relationship between the median R 2 m and the MAF value was then plotted to characterize the relationship.
The impact of the number of nonmissing data points at a marker on the marker's imputation accuracy was assessed for each datasetimputation method combination using data from all 10 replicates and versions NA20, NA50 and NA70 combined. For each marker, the number of nonmissing data points was rounded to the nearest factor of 5, and for each value the median R 2 m was calculated. To determine the impact of the LD level with other markers on the imputation accuracy, markers were first classified as markers in low LD with all other markers or markers in at least moderate LD with at least one other marker. Markers whose highest r 2 statistic was less than 0.5 were considered to be in low LD with all other markers. A marker that had at least one r 2 statistic greater than or equal to 0.5 was considered to be in at least moderate LD with at least one other marker. The median R 2 m of markers in low LD and of markers in at least moderate LD with at least one other marker was calculated. The ratio of R 2 m for markers in low LD to the R 2 m for markers in at least moderate LD was then examined.
To assess the effect of the genetic distance between an individual and its closest relative on the individual genotype imputation accuracy, the Euclidian distance was calculated for each pair of individuals and the R 2 i of each dataset was measured for each simulated dataset and imputation method combination. The mean R 2 i values across all replicates, R 2 l , were averaged across versions NA20, NA50, and NA70 of a given dataset-imputation method combination to calculate an overall mean R 2 i for each individual which is referred to as R 2 l . The Euclidian distance between each individual and its closest relative, rounded to the nearest whole number was plotted against the median R 2 l to examine the relationship.
The relationship between PEV for the genetic values and the R 2 l was also examined. An individual's PEV is a measure of genetic connectedness to the other individuals (Kennedy and Trus 1993) where an individual's connectedness is determined by the number and strength of the genetic relationships between that individual and the other individuals in the dataset. For example, a low PEV indicates high connectedness and high degree of genetic relationship.
To measure an overall PEV value for each individual, a vector of PEVs was calculated for each marker using the mixed model equations (Searle et al. 1992) implemented in the R package "rrBLUP" (Endelman 2011). The genetic and error variance components were estimated using restricted maximum likelihood, and the genomic relationship matrix, excluding the response variable marker, was used as the covariance matrix between genotypes. The sum of the PEV vectors across all markers was used as the overall PEV vector. Because PEV is a reflection of the number and strength of the genetic relationships between individuals, it is expected to be a useful indicator for how well and individual's missing data can be imputed using all other individuals as a reference.

GS accuracy calculation
All 10 simulations of missing data versions NA20, NA50, and NA70 of the WW, SW, DTM, NAB, and SRRW marker sets were imputed with each of the imputation methods: MNI, kNNI, SVDI, EMI, and RFI. Then, each of the 10 replicates of the marker set-imputation method combinations was used to calculate the 10-fold cross validation GS accuracy for both Ridge-Regression (Whittaker et al. 2000) and Bayesian LASSO (de los Campos et al. 2009), see File S1. GS accuracies are computed as the Pearson's correlation between the phenotype estimated breeding values and the genomic-estimated breeding values. The mean accuracy for each marker set-imputation method-prediction model combination was computed. GS accuracies were also computed using version NA0 of the WW, SW, DTM, and NAB genotypic data.

LD between markers
For each dataset, the LD between marker pairs was quantified using the r 2 statistic. Markers that had at least one other marker associated with r 2 $ 0.5 were considered to be in at least moderate LD with at least one other marker. In the WW, SW, DTM, and NAB datasets, 62%, 74%, 12%, and 69% of the markers had at least one other marker in at least moderate LD, respectively. Comparatively, LD between markers was high in the SW, NAB, and WW datasets and much lower in the DTM dataset.

Imputation accuracy
The imputation accuracy reported as the median R 2 m is shown in Figure 1 for kNNI, SVDI, RFI, and EMI. For all dataset-imputation method combinations, R 2 m values were non-normal, and there were  Table 2. The population with the highest median R 2 m for each of the levels of missing data were the NAB population, whereas the lowest imputation accuracies were observed with the DTM population. As expected, median R 2 m values always decreased as the level of missing data increased. RFI always produced the highest accuracies; kNNI generally produced the second highest accuracies, followed by EMI and SVDI. The rankings were slightly different for the DTM dataset, where RFI was most accurate followed by EMI, SVDI, and kNNI. The rankings of the methods for each dataset according to the median percent correct are the same as those according to the median R 2 m ; however, the median percent correct values could not be compared across datasets because percent correct values are influenced by the MAF which differs among datasets.

Computational time
Large differences in the computational requirements for the imputation methods were observed (Table 3). kNNI, SVDI, and EMI required relatively little computation time on average, while RFI required at least 95x, 760x, and 65x more computation time than kNNI, SVDI, and EMI, respectively. For SVDI and kNNI, the computation time required for determining optimal k values was not included in the estimates of the average computational time because the computation time for optimal k estimation depends on the method used for estimation. The 10-fold cross validation approach that we used to estimate optimal k values for SVDI and kNNI requires approximately 50 runs of the SVDI and kNNI respectively. If the time required to estimate optimal k values for SVDI and kNNI was included in the total computational time, EMI would be the fastest of the four imputation methods.

Factors affecting imputation accuracy
MAF: For all datasets, R 2 m values for markers with MAF , 0.1 were low compared with that of markers with MAF . 0.1; however, the relationship between MAF and R 2 m for markers with MAF . 0.1 was different for each dataset (Figure 2). In general, R 2 m increased as MAF increased as long as MAF , 0.4; however, with the NAB dataset ( Figure  2D) there was no relationship between MAF and R 2 m for MAF . 0.1. Accuracy in terms of percent correct had a strong negative linear relationship with the MAF across all imputation methods and datasets. Markers with lower MAF values tended to have higher percent correct values (data not shown).

Number of nonmissing data points:
With almost all datasetimputation method combinations, as the number of non-missing data points increased, the R 2 m levels increased in a linear fashion ( Figure 3). The strength of this linear relationship was similar for all imputation methods; however, with the DTM dataset, R 2 m for kNNI and SVDI were close to zero regardless of the number of nonmissing data points.

LD between markers:
The ratio of the median imputation R 2 m for markers with no other markers in moderate LD to the median imputation R 2 m for markers with at least one other marker in moderate LD was always less than one (Table 4), indicating that the imputation accuracy for markers without markers in moderate LD was always lower than that for markers that had at least one other marker in moderate LD. Across all datasets, the R 2 m ratios for the two classes of markers was much smaller for kNNI compared to the other imputation methods, indicating that the imputation accuracy of kNNI was more strongly influenced by the level of LD between markers compared to the other methods. With the WW, SW, and NAB datasets the R 2 m ratios for the two classes of markers was similar for SVDI, RFI, and EMI indicating that the accuracy of these three methods is influenced by the level of LD between markers to a similar degree. However, with the DTM dataset, the R 2 m ratio for the two classes of markers was closer to one for SVDI compared to the other methods, indicating that for this dataset, the accuracy with SVDI was less affected by the LD between markers, compared to the other methods.
n Distance from the closest relative and PEV: Regardless of the dataset or the imputation method, the smaller the distance between an individual and its closest relative, the higher the R 2 l (Figure 4). One exception was observed with the DTM dataset, where for kNNI there was no relationship between the distance between an individual and its closest relative and R 2 l . We observed very similar trends between R 2 l and the overall PEV ( Figure S2). As an individual's PEV increased, indicating a decrease in the strength and number of genetic relationships between that individual and all other individuals, its R 2 l decreased in all cases except when the DTM dataset was imputed with kNNI.

Effect of imputation method on GS accuracy
In nearly all cases, GS accuracies did not differ greatly from one imputation method to another, with the exception of MNI, which sometimes led to much lower accuracies compared to all other methods when the NA70 dataset version was used ( Figure 5 and 6). Overall, GS accuracies were least affected by the imputation method for dataset version NA20, and most affected by the imputation method for dataset version NA70. The relative performance of each method in terms of GS accuracy after imputation depended on the dataset, and dataset version; however, RFI consistently performed well across all datasets. For the WW datasets, the relative performance of the imputation methods in terms of GS accuracy was inconsistent across the four traits tested.
For, a given dataset and dataset version, the rank of each method based on R 2 m , was not consistent with the rank based on GS accuracy using RR or BL post-imputation. The rank of the imputation methods, however, was consistent between the two different GS models. We also found that including rather than removing 'sparse' markers, those with large amounts of missing data, nearly always led to higher GS accuracies (methods and results described in File S1, Table S2, and Figure S3), especially when RFI, kNNI, or EMI were the imputation methods used ( Figure S4).

Imputation accuracy
This study found that map-independent imputation methods other than MNI can be surprisingly accurate, especially when LD between markers is high and the genotyped individuals are related. RFI was the most promising method overall because of its consistently high performance in terms of imputation accuracy and subsequent GS accuracy; however, it was the most computationally intensive method evaluated. kNNI, although less accurate than RFI, may be a good alternative to RFI if there are computational limitations to completing the imputation. It is likely that RFI and kNNI produced comparable levels of accuracy because both use a similar model free approach for imputation that involves weighting a selected set of k important variables according to a distance metric (Lin and Jeon 2006). The weighted average of these variables is the predicted value of the variable of interest. For kNNI the distance metric was the Euclidean distance and k was a fixed number across all variables. For RFI, the k important variables and their weights are determined by the splitting scheme of the tree that is determined using the response variable. The The median R 2 m obtained for a given MAF rounded to the nearest tenth is plotted for each dataset: (A) WW; (B) SW; (C) DTM; and (D) NAB. Each color and symbol represents a different imputation method: kNNI, orange triangles; SVDI, purple squares; RFI, red circles; and EMI, blue crosses.
increased accuracy but greater computational burden of the RFI method compared with kNNI is due to its adaptive weighting of variables that takes into account the response variable.
A possible reason that EMI and SVDI were less accurate than RFI and kNNI is that the genotypic datasets that we used may have violated multivariate normality, an underlying assumption for EMI and SVDI. Alternatively, EMI and SVDI may not have been as effective at ignoring uninformative predictors. If true, linear regression based imputation methods involving variable selection could be as accurate as kNNI or RFI. However, due to multicollinearity, attempts to test imputation based on subset selection methods such as stepwise regression were not successful. Regression imputation using variable selection methods which can cope with multicollinearity, such as the least absolute shrinkage and selection operator (Tibshirani 1996), would be interesting to test in future studies.
EMI performed consistently better than SVDI which is likely because EMI incorporates all the marker data as predictors whereas SVDI first used a data reduction step, potentially eliminating useful information. SVDI may have outperformed EMI if the datasets had a greater rate of genotyping error because it is expected to better cope with noisy data (Troyanskaya et al. 2001).
For all methods, average median imputation accuracies on an individual genotype basis R 2 l were not always homogenous across population sub-groups as illustrated in Figure S5, which shows individuals plotted according to the first two principal components of their marker genotypes and color coded according to their n Table 4 Ratios of median R 2 m of markers having no markers in moderate linkage disequilibrium (LD) to that of markers with at least one other marker in moderate LD imputation accuracy. With the DTM and WW datasets, small subgroups of individuals that clustered together according to the first two principal components of marker genotypes tended to have similar ranges of accuracy. However, with the SW and NAB datasets R 2 l was relatively homogenous across population subgroups. An association between R 2 l and population subgroup is undesirable because it may create or worsen an association between GS accuracy and population sub-group. Using large datasets with minimal population structure for imputation and GS is advocated to avoid heterogeneity of imputation and GS accuracies across subgroups of individuals.
Population structure may also lead to increased imputation accuracy for markers with high levels of population subdivision (Iwata and Jannink 2010) because an individual's allelic state can be predicted largely based population sub-group alone. Accuracy levels for datasets with many markers highly subdivided by population may be high largely because of structure; we therefore calculated R 2 m excluding markers with high levels of population subdivision as indicated by their Fst values, where high Fst indicates high population subdivision (for methods, see File S1). For markers with MAF . 0.1, on average, R 2 m excluding markers with the 25% highest Fst values were 0.9, 1.17, 1.02, and 0.9 times those of overall R 2 m for the WW, SW, DTM, and NAB datasets, respectively. Thus, for the WW and NAB datasets, the high imputation accuracies we observed may have been in small part due to population structure.
Comparing our imputation accuracy results with those of other studies is difficult because each study uses different populations of different sizes, levels of missing data, MAF distributions, and levels of LD between markers. In addition, accuracy reported as percent correct cannot be compared across datasets with different MAF distributions. Nevertheless, we assume that map-dependent imputation methods would outperform the map-independent methods that we evaluated (given the availability of an accurate genetic or physical map) because physically linked markers are used to predict missing values. These physically linked markers should be more reliable predictors compared to markers that are in LD but may not be physically linked. As genetic and physical maps develop for wheat and barley the assumption that map-dependent methods would outperform the map-independent methods can be tested.

Factors affecting imputation accuracy
Markers with very low MAF had low R 2 m values. There are two possible explanations for this observation. First, because of the way R 2 m is calculated, a single imputation error has a much larger negative impact on the R 2 m for markers with lower MAF values ( Figure S6). Thus, it is harder to achieve high R 2 m for markers with a low MAF. Second, individuals with the minor allele at a given marker are not well represented, making their marker genotype more difficult to predict. A similar relationship between MAF and R 2 m was also found by studies by Iwata and Jannink (2010) and Li et al. (2011), which used mapdependent imputation methods. Unlike R 2 m , imputation accuracy in terms of percent correct had a negative linear relationship with MAF (data not shown), this is because markers with lower MAF can always be imputed with a reasonably high percent correct based on the marker mean alone. Other studies of map-dependent imputation methods report a negative relationship between MAF and percent correct (Pei et al. 2008;Hickey et al. 2012). The number of nonmissing data points, analogous to reference panel size in other studies was found to positively impact the R 2 m . This finding is consistent with other studies in which researchers tested the effect of reference panel size on the imputation accuracy using mapdependent methods (Pei et al. 2008;Druet et al. 2010;Li et al. 2010). For RFI, EMI, and SVDI, which involve a model training step, fewer missing data points means that more individuals are available for model training. With kNNI, a smaller number of nonmissing data points at a given marker leads to a more accurate estimate of its distance from all other markers. However, with the DTM set there was no trend between accuracy and the number of nonmissing data points with kNNI. This was because accuracy with kNNI for this dataset was very low overall.
The presence of one or more markers in moderate LD (r 2 statistic $ 0.5) was a more important factor for kNNI compared with RFI, EMI, and SVDI because kNNI bases its predictions on a fixed number of close markers, whereas RFI, EMI, and SVDI use information from all markers in the dataset to generate predicted values for the missing data points. The LD between markers on a whole dataset basis also appeared to be an important factor affecting the R 2 m of all methods because accuracies with the DTM dataset, which had low levels of LD between markers overall, were much lower than accuracies with the WW, SW, and NAB datasets. Other publications that have evaluated the effect of LD on imputation accuracy for map-dependent methods have found similar trends (Pei et al. 2008;Hickey et al. 2012).
We found that imputation accuracy on an individual genotype level was negatively correlated with the distance from the closest relative in the dataset, and the PEV, which is an indication of the relationship between an individual and other genotypes. A similar relationship between imputation accuracy and relationship has been found by other studies of map-dependent imputation methods Zhang and Druet 2010;Hickey et al. 2012). It is clear that to ensure effective imputation, the dataset to be imputed should contain related individuals. If the dataset is suited for GS, it is likely that the individuals are already related. However, to increase the chances Figure 5 GS accuracy obtained using ridge regression after imputation. Mean GS accuracies obtained using the dataset versions NA0, NA20, NA50, having up to 0%, 20%, 50%, and 70% missing data per marker, respectively, imputed with either MNI (black stars), kNNI (orange triangles), SVDI (purple squares), EMI (blue crosses) and RFI (red circles) are shown for (A) WW-yield, (B) WW-height, (C) WW-protein, (D) WW-days to heading, (E) DTM, (F) SW, (G) NAB, and (H) SRRW datasets. Each plot has a different y-axis range. Error bars depict SE.
that an individual will have close relatives in the dataset, all available genotypic data for the germplasm pool of interest should be combined before imputation.

GS accuracy
The GS accuracies that we observed may be sufficiently high to lead to increased rates of genetic gain compared to phenotypic selection (PS), depending on the accuracy of PS and the selection cycle duration of both PS and GS. It is important to note that all GS accuracies reported for a given dataset are global estimates across all potential subpopulations. Based on other studies evaluating GS accuracies within and across subpopulations (Zhao et al. 2011, Heslot et al. 2012, Windhausen et al. 2012, this global accuracy estimate may be greater than the accuracy measured within individual subpopulations. Effect of imputation method on the GS accuracy Improved accuracy of GS after application of map independent imputation methods was another important finding of this study. Based on our results, unordered markers with missing data can be included in the dataset to improve accuracy through imputation with RFI, kNNI, EMI, or even SVDI rather than MNI. However, for datasets with low levels of missing data (up to 20% per marker), imputing with MNI is sufficient. Although our results do not support removing markers with high levels of missing data prior to GS, in many datasets markers with low levels of missing data may be sufficient to saturate the genome. With the datasets used in this study, the average number of markers with up to 20% and 50% missing data were 18 to 37 and 99 to 186 respectively, and these reduced marker sets were not sufficient to saturate the genome. Thus, including markers with larger amounts of missing data led to improved GS accuracies. Interestingly, a low median R 2 m was not reflective of the merit of imputation prior to GS. The median R 2 m for the datasets with up to 70% missing data per marker were the lowest of all the missing data levels; however we saw the greatest gain in GS accuracy from kNNI, SVDI, EMI, or RFI relative to MNI with this level of missing data. This was especially apparent for the DTM dataset, which had a median R 2 m near zero for most methods when there was up to 70% missing data per marker. However, RFI on this dataset produced GS model accuracies 1.3 times Figure 6 GS accuracy obtained using Bayesian least absolute shrinkage and selection operator after imputation. Mean GS accuracies obtained using the dataset versions NA0, NA20, and NA50, having up to 0%, 20%, 50%, and 70% missing data per marker respectively, imputed with either MNI (black stars), kNNI (orange triangles), SVDI (purple squares), EMI (blue crosses) and RFI (red circles) are shown for (A) WW-yield, (B) WW-height, (C) WW-protein, (D) WW-days to heading, (E) DTM, (F) SW, (G) NAB, and (H) SRRW datasets. Each plot has a different y-axis range. Error bars depict SE. greater than those achieved when MNI was used before GS. Surprisingly, the most accurate imputation method was not always the method that gave the greatest GS accuracy. This may be caused by nonrandom imputation errors. If some imputation errors are similar for related individuals, these nonrandom errors may able to capture some genetic relationships in the GS model. The idea that the imputation errors may capture some genetic relationships was suggested by a study by Weigel et al. (2010).
This study has important implications for species that lack a reference genome, complete reference map, and predesigned high-throughput genotyping platforms. First, unordered markers can be imputed with high levels of accuracy, and even greater accuracies may result if additional reference genotypes can be added to the dataset prior to imputation. Based on the results of this study, if a large number of marker genotypes are produced (so that markers are in LD with each other), and the population contains individuals with some genetic relationship, missing data can be imputed with reasonable accuracy even if the level of missing data are high; up to 70%. Future work to improve upon and reduce the computational burden of the most promising methods in this study, RFI and kNNI, would be especially useful if these methods are to be used widely. The second implication of this study is that a large proportion of missing data in dense marker sets is not a major concern for GS. As long as the marker density is sufficiently high, the accuracy does not appear to be strongly negatively affected. In cases where missing data does negatively impact the GS accuracy imputation using a method other than MNI before GS model training and validation can help improve the accuracy. Overall, mapindependent imputation shows promise for the feasibility of applying GS, enabled by emergent sequence-based genotyping technologies, to almost any species regardless of the availability of pre-existing genotyping resources.


extfieldsofstudy: Chemistry Medicine

text:
Acarbose Reduces Blood Glucose by Activating miR-10a-5p and miR-664 in Diabetic Rats

MicroRNAs (miRNAs) are non-coding RNA molecules involved in the post-transcriptional regulation of a large number of genes, including those involved in glucose metabolism. Acarbose is an α-glucosidase inhibitor that improves glycemic control by decreasing the intestinal absorption of glucose, thereby decreasing the elevation of postprandial blood glucose. However, acarbose is poorly absorbed into the blood stream from the gut. Therefore, the exact mechanisms by which acarbose affects glucose metabolism are unclear. This study investigated the effect of acarbose on glucose metabolism in diabetic rats and tested the hypothesis that acarbose acts directly through miRNA-regulated expression in the intestinal epithelium. Rats were divided into four groups: a control group, a diabetic group (DM), a low dose of acarbose group (AcarL) and a high dose of acarbose group (AcarH). Ileum samples were analyzed using miRCURY LNA™ microRNA Array, qPCR and immunohistochemistry. We found that 8-week treatment with acarbose significantly decreased fasting blood glucose. Oral glucose tolerance tests (OGTT) showed that blood glucose was significantly reduced in the AcarL and AcarH groups at 30 min, 60 min and 120 min after oral glucose administration. We found that miR-151*, miR-10a-5p, miR-205, miR-17-5p, miR-145 and miR-664 were up-regulated in the AcarH group, while miR-541 and miR-135b were down-regulated. Through target gene analysis, real time PCR and immunohistochemistry verification, we found that these miRNAs suppressed the expression of proinflammatory cytokines [IL6 (interleukin 6) and TNF (tumor necrosis factor)] and mitogen activated protein kinase 1 (MAPK1). Our data suggest that acarbose can improve blood glucose in diabetic rats through the MAPK pathway and can down-regulate proinflammatory factors by activating miR-10a-5p and miR-664 in the ileum.

Introduction
Diabetes mellitus is one of the most common chronic diseases worldwide, and continues to increase in incidence and significance, as changing lifestyles lead to reduced physical activity, and increased obesity. Type 2 diabetes mellitus is an emerging worldwide health problem, with the number of global cases of type 2 diabetes projected to double to 350 million by the year 2030 [1]. Diabetes is an independent risk factor for cardiovascular disease [2,3] and is the leading cause of morbidity and mortality in the developed world [4][5][6].
Acarbose is an a-glucosidase inhibitor that delays the digestion of complex carbohydrates and disaccharides to absorbable monosaccharides by reversibly inhibiting a-glucosidases within the intestinal brush border, thereby attenuating postprandial blood glucose peaks [7]. Clinical trials have demonstrated that acarbose generally improves glycemic control in patients with diabetes mellitus that can be managed by diet alone, or in combination with other antidiabetic therapies, as evidenced by decreased postprandial plasma glucose and glycosylated hemoglobin. It does not appear to directly alter insulin resistance, but it may lower postprandial plasma insulin levels. However, the bioavailability of acarbose is low [8], which is attributed to its poor aqueous solubility.
MicroRNAs (miRNAs) are short (21-23 nucleotides), endogenous, non-coding RNA molecules. miRNAs regulate gene expression by imperfect base pairing with the 39-untranslated regions of mRNAs, resulting in mRNA decay or translational repression [9]. miRNAs have distinct spatial and temporal expression patterns in cells and tissues and regulate several processes, including hematopoiesis, development, cell differentiation, proliferation and apoptosis [10,11]. They are implicated in several diseases, including diabetes.
We therefore hypothesized that acarbose directly alters the intestinal expression of miRNAs to regulate glucose metabolism. To provide molecular evidence for this mechanism, we used a rat model of type 2 diabetes to investigate differential miRNA expression in rat intestines after treatment with acarbose. and then were administered a single dose of streptozotocin (STZ, 50 mg/kg, tail vein) formulated in 0.1 mmol/L citrate buffer, pH 4.5 (Sigma-Aldrich, MO, USA). One week after the STZ injection, the random blood glucose level of the diabetic rats was measured to confirm hyperglycemia. Random blood glucose above 16.7 mmol/L was used to define rats as diabetic. Diabetic rats were fed a high-fat diet throughout the experiment. Diabetic rats with a similar degree of hyperglycemia were randomly divided into three groups: vehicle, low dose acarbose (AcarL), and high dose acarbose (AcarH) groups (n = 10, in each group). The typical human daily dose of acarbose is 300 mg/60 kg body weight. According to the formula: d rat = d human 6 0.71/0.11 [13], the corresponding dose of acarbose for rats is 32.28 mg/kg per day. Therefore, we selected 30 and 60 mg/kg per day as low and high dosages, respectively. The control (n = 10) and the diabetic group received 0.5% saline, whereas the AcarL and AcarH groups were given acarbose (Bayer Health Care Co., Germany) at doses of 30 and 60 mg/kg in a 0.5% saline solution, respectively. The drug was administered once daily for 8 weeks using a gastric gavage. All animals were housed in an environmentally controlled room at 25uC with a 12 h light-dark cycles and were given free access to food and water throughout the experimental period. Fasting animals were allowed free access to water. After 6 weeks of treatment, an oral glucose tolerance test (OGTT) was performed. After 8 weeks of treatment, blood samples were taken from rats after anesthesia. The rats were then sacrificed. Some terminal ileums were collected for performing the microarray and quantitative real-time reverse transcription PCR (qRT-PCR) analysis. Other terminal ileums were fixed in 10% neutralized formalin for immunohistochemical staining. All procedures involving animals were approved by the animal care and use committee of the Peking Union Medical College Hospital (Beijing, China, MC-07-6004) and were conducted in compliance with the Guide for the Care and Use of Laboratory Animals, 8th ed., 2011. All surgeries were performed under sodium pentobarbital anesthesia, and all efforts were made to minimize suffering.

Measurement of Body Weight and Fasting Blood Glucose
Body weight was measured every 2 weeks. The 6-h fasting blood glucose (FBG) level was measured every 2 weeks using a Contour TS glucometer (Bayer) with blood from a tail bleed.

Serum IL6 and TNF-a Level Analysis
At week 8, blood samples were collected after euthanasia and centrifuged at 1000 g for 10 min. Serum was stored in aliquots at 280uC to assay serum interleukin 6 (IL6) and tumor necrosis factor a(TNF-a). Serum IL6 and TNF-a levels were measured by enzyme-linked immunosorbent assay (ELISA, Abcam, UK).

miRCURY LNA TM microRNA Array Experiment
The miRCURY LNA TM miRNA array system contains 3100 capture probes, covering all rat microRNAs (388 miRNAs) that have been annotated in miRBase 18.0 as well as all viral microRNAs associated with rats. Total RNA from the iliem of AcarH group and DM group was harvested using TRIzol (Invitrogen) and an miRNeasy Mini Kit (QIAGEN) according to the manufacturers' instructions. After the RNA was quantified using a NanoDrop 1000, the samples were labeled using a miRCURY TM Hy3 TM /Hy5 TM Power labeling kit and hybridized on a miRCURY TM LNA Array v.18.0 (Exiqon). Following washing, the slides were scanned using an Axon GenePix 4000B microarray scanner.

Gene Array Data Analysis
Normalization was performed with a per-chip 50 th percentile method that normalizes each chip on its median, allowing for comparison among chips.

miRNA Quantitative Real-time PCR (qRT-PCR)
Total RNA (5 ng) was reverse-transcribed using a Taqman TM MicroRNA Reverse Transcription kit (Applied Biosystems) and the miRNA-specific reverse-transcription primers provided with TaqMan TM MicroRNA Assays (Applied Biosystems). For reversetranscription, a PTC-225 Peltier Thermal Cycler (MJ Research Inc., Waltham, Massachusetts) was used. The reaction conditions were 16uC for 30 min, 42uC for 30 min and 85uC for 5 min. The generated miRNA-specific cDNA was amplified using a Taq-Man TM Universal PCR master mix II (Applied Biosystems) and the respective specific probe provided with TaqMan TM Small RNA Assays (Applied Biosystems). PCR was performed using a CFX-96TOUCH (Bio-Rad) detection system. Amplification was performed at 95uC for 10 min, followed by 40 cycles of 95uC for 15 sec and 60uC for 60 sec. U6 small nuclear RNA was used as an endogenous control. The fold change in miRNA level was

Target Gene Validation by qRT-PCR
For the validation of miRNA target genes, qRT-PCR analyses of RNAs were performed using SYBR Green. Each qRT-PCR assay was repeated using three biological replicates, and each analysis consisted of three technical replicates. Before PCR analysis, each sample of total RNA was treated with RNase-free DNase (Qiagen, Valencia, CA, USA). RNA was reversetranscribed by Superscript II (Invitrogen, CA, USA). The primers were designed using the Applied Biosystems (Foster City, CA, USA) Primer Express TM design software. Primers were purchased from Applied Biosystems (Table 1). Using the ABI Prism 7700 Sequence Detection System, the following reaction conditions were used: an initial denaturation at 48uC for 30 min, followed by 95uC for 15 min, and then 40 cycles of 95uC for 15 sec, and 55uC for 1 min, and a final unlimited 4uC hold. The sequences of the primers are listed in Table 1. The signal of the housekeeping gene Gapdh (glyceraldehyde-3-phosphate dehydrogenase) was used for normalization. The relative quantification of the mRNA between the acarbose treated groups and DM rats was calculated using the comparative Ct method.

Immunohistochemical Staining
Ileum from the AcarH and DM groups (n = 6 in each group) were fixed in 10% neutral buffered formalin, cast in paraffin, sliced into 4-mm sections and placed onto microscope slides. After the removal of the paraffin by xylene and dehydration by graded alcohol, the slides were immersed into distilled water. Ileum sections were then transferred into a 10-mmol/L citrate buffer solution (pH 6.0) and heated at 80uC for 5 min for antigen retrieval. After washing, 3.0% peroxide was applied for 20 min to block the activity of endogenous peroxidase. To avoid nonspecific staining, the sections were incubated in blocking solution (5% BSA) for 1 h at room temperature, followed by treatment with mouse monoclonal anti-TNF-a antibody (1:50, Abcam Inc., UK), rabbit polyclonal anti-IL6 antibody (1:100, Abcam lnc., UK), or rabbit polyclonal anti-MAPK1 antibody (1:100, Abcam Inc., UK) where indicated, overnight at 4uC. Negative control sections were stained under identical conditions by substituting the primary antibody with equivalent concentrations of normal rabbit IgG. After washing with phosphate-buffered saline, the slides were incubated with a labeled streptavidin biotin reagent, following the manufacturer's instructions. Immunoreactive products were visualized with the DAB reaction. Sections were counterstained with hematoxylin for 15 sec. Brownish yellow granular or linear deposits were interpreted as positive areas. Three observers blinded to the clinical information evaluated the immunohistochemical staining scores independently. Staining intensity was graded semi-quantitatively using the H-SCORE [15] which was calculated using the following equation: H-SCORE = gPi (i +1), where i is the intensity of staining with a value of 1, 2 or 3 (mild, moderate, or strong, respectively) and Pi is the percentage of epithelial cells stained with different intensity, varying from 0% to 100%. The results are expressed as the mean 6 SE. Differences between the groups were statistically analyzed with a one-way analsis of variance (ANOVA). A P value of ,0.05 was considered significant.

Statistical Analysis
All results are expressed as the mean 6 standard deviation (SD). Statistical analyses were performed with ANOVA followed by Student's t test. P,0.05 was considered statistically significant. Analyses were performed with SPSS 11.0 (SPSS, Inc., Chicago, IL, USA).

Acarbose had No Effect on the Body Weight of DM Rats
The mean body weight of diabetic rats significantly decreased compared to the control rats at week 2 (P,0.05), week 4 (P,0.01), week 6 (P,0.01) and week 8 (P,0.01). No significant differences were noted between the DM group and the acarbose-treated groups ( Figure 1A).

Acarbose Regulated the Glucose Tolerance of DM Rats
The blood glucose levels of the DM group were higher than those of the control group before and at 30 min (P,0.01), 60 min (P,0.01) and 120 min (P,0.01) after oral glucose administration. Blood glucose levels of the acarbose-treated groups significantly decreased at 30 min, 60 min and 120 min after oral glucose administration (P,0.05, Figure 2).

Acarbose Reduced Serum IL6 and TNF-a in DM Rats
To determine whether acarbose suppressed inflammatory mediators in DM rats, the levels of serum IL6 and TNF-a were determined. The levels of serum IL6 and TNF-a were significantly elevated in DM rats (P,0.01). Acarbose significantly suppressed serum IL6 (P,0.05) and TNF-a (P,0.01) after the 8-week treatment (Figure 3).

miRNA qRT-PCR
To validate the microarray results, all of the miRNAs were selected for qRT-PCR quantification. All of the miRNA expression levels obtained by qRT-PCR were similar to those observed by microarray analysis (Figure 5).

Identification of the Target Genes of the Differentially Expressed miRNAs Using Bioinformatics Analysis
We next performed bioinformatics analysis to identify the target genes of the differentially expressed miRNAs. Together, the eight miRNAs had 189 validated target genes in the miRWalk database (Table 3). By using the miRWalk and KEGG databases, these 189 validated genes were implicated in the TGF-b signaling pathway, the MAPK signaling pathway, the Wnt signaling pathway and to be involved in the cytokine-cytokine receptor interaction (Table 4).

mRNA qRT-PCR
To determine whether local gene expression of inflammatory markers in the ileum were altered by acarbose treatment, the mRNA expression of IL6, Mapk1 and TNF were determined by qRT-PCR. The expression of Il6, Mapk1 and Tnf was downregulated significantly in the AcarH group (P,0.05, Table 5).

Immunohistochemical Staining
To confirm that the protein expression of inflammatory markers were altered in the ileum of DM rats treated with acarbose, immunohistochemistry analyses for TNF-a, IL6 and MAPK1 were performed on ileum tissues. In the AcarH group, there was a statistically significant decrease in the immunoreactivities of TNFa, IL-6 and MAPK1 ( Figure 6).

Discussion
In this study, the administration of acarbose to diabetic rats significantly reduced blood glucose and altered the expression of miRNAs in the intestine. The miRNAs altered were found to target genes in inflammatory pathways, including IL6, Mapk1 and Tnf. In addition, serum IL-6 and TNF-a levels were shown to be significantly decreased in the AcarH group. These results indicate that acarbose reduces inflammation, potentially via miRNAregulated signals in the intestine. To the best of our knowledge, this is the first report of glucose moderation of acarbose by miRNAs through the inflammatory pathway.
Acarbose is poorly absorbed into the blood stream from the gut. Therefore, we hypothesized that acarbose has a direct effect on the gut to achieve glucose moderation and inflammation reduction. To investigate the effect of acarbose on miRNA expression, potential miRNA targets in the terminal ileum were identified. In our miRNA gene array analysis, eight miRNAs showed significant changes in the AcarH group compared to the DM group. The other miRNAs had 189 target genes that have been experimentally verified. A bioinformatic analysis of all differentially expressed miRNAs in the ileum of the AcarH group, suggested that these genes were involved in the TGF-b signaling pathway, the MAPK signaling pathway, and the Wnt signaling pathway as well as the cytokine-cytokine receptor interaction. qRT-PCR and immunohistochemical staining verified that the expression of Il6, Tnf and Mapk1 were down-regulated in the AcarH group. Il6 and Tnf were identified as target genes of miR-10a-5p, and Mapk1 was one of the target genes of miR-664. These findings suggest that acarbose can reduce the expression of inflammatory cytokines in the ileum through miR-10a-5p and miR-664.
Individuals with type 2 diabetes display the features of lowgrade inflammation [16]. Cytokines are a group of pharmacologically active low molecular weight proteins possessing autocrine and paracrine effects [17]. Moreover, they are products and effectors of the inflammatory and immune systems [18]. Chronic hyperglycemia modulates immune responses and triggers inflammation through the activation of Toll-like receptors (TLRs), leading to an increase in the proinflammatory cytokines, IL1, IL6, TNFa and interferon c(IFNc), and a decrease in the levels of interleukin 10 (IL10). Proinflammatory cytokines inhibit insulin signals by activating receptor kinase inhibitors, stimulating NFkB to induce pancreatic b cell damage and apoptosis. IL6 can affect glucose homeostasis and metabolism through multiple pathways [19]. IL6 can also increase plasma plasminogen activator inhibiotor (PAI-1) and c-reactive protein (CRP) levels [20] and inhibit glucose-stimulated insulin secretion in islets [21]. Acarbose may upregulate miR-10a-5p to suppress the expression of IL6 in the gut, thereby moderating glucose homeostasis and inflammation in diabetic rats.  TNF-a is another target gene of miR-10a-5p. TNF-a reduces glucose transporter 4 in cells and decreases glucose uptake. TNF-a also impairs insulin action interfering with insulin signaling by binding to TNF-a receptors on muscle cells or hepatocytes [22]. Administering soluble TNF-a receptors in obese rats neutralizes TNF-a and improves insulin sensitivity [23]. Continuous exposure of 3T3-L1 adipocytes to TNF results in a decrease in GLUT4 and insulin receptor (IR) gene expression [24]. Thus, acarbose may moderate glucose metabolism and proinflammation by upregulating miR-10a-5p and decreasing the expression of TNF-a.
Using the KEGG analysis, we observed that the primary pathway through which miRNAs targeted genes in the acarbose  group was the MAPK signaling pathway. qRT-PCR verified this result. MAPK1 was identified as a target gene of miR-664, and TNF-a is interconnected with MAPK pathways. Waetzig et al. found that SB 203580, a p38 inhibitor, significantly reduces the secretion of TNF-a in the mucosa [25]. Hollenbach et al. found the SB203580 reduces the mRNA levels of proinflammatory cytokines (e.g.,TNF-a, IL2, IL18) in the gut of BALB/c mice [26]. Thus, we propose that acarbose increases the expression of miR-664, leading to decrease Mapk1 expression, which ultimately reduces TNF-a in the gut.
In conclusion, acarbose can regulate glucose metabolism through the MAPK pathway and can suppress of proinflammatory cytokines by increasing miR-10a-5p and miR-664 in the ileum. These results provide molecular information for the future investigation of the mechanisms by which acarbose regulates glucose metabolism in DM rats at the miRNA level. Furthermore, these results could be important in devising mechanism-based and targeted therapeutic strategies for DM. However, future studies using human samples must be carried out to validate our findings and to further elucidate the miRNAs regulatory networks that are affected by acarbose in the treatment of type 2 diabetes.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Psychology

text:
A mobile genetic element with unknown function found in distantly related viruses

Background The genetic element s2m seems to represent one of very few examples of mobile genetic elements in viruses. The function remains obscure and a scattered taxonomical distribution has been reported by numerous groups. Methods We have searched GenBank in order to identify all viral accessions that have s2m(−like) sequence motifs. Rigorous phylogenetic analyses and constrained tree topology testing were also performed in order to investigate the apparently mobile nature of s2m. Results The stem-loop s2m structure can be found in four families of + ssRNA viruses; Astroviridae, Caliciviridae, Picornaviridae and Coronaviridae. In all of these virus families, with the possible exception of Caliciviridae, multiple gains and/or losses of s2m would have to be postulated in order to explain the distribution of this character. Conclusions s2m appears to be a mobile genetic element with a unique evolutionary history in all of the four virus families where it can be found. Based on our findings and a review of the current literature on s2m, a hypothesis implying an RNAi-like function for the s2m element is also outlined.

Background
In 1997, a 43 basepair, conserved sequence motif was described near the 3' end of members of the Astroviridae family [1]. The genetic element corresponded to the second most 3' stem-loop structure (stem-loop II) in human astroviruses and was subsequently named s2m [2]. The sequence motif has later been found in three other virus families; Caliciviridae, Picornaviridae and Coronaviridae [3]. The distribution of s2m seems to be limited to positive-sense, single-stranded RNA (+ssRNA) viruses and the element is always located near the 3' end of the genome. Most commonly, s2m can be found downstream of the last reading frame, but there are also instances where it appears that the final stop codon is part of the motif itself [2]. Recently, several examples of viruses containing two copies of s2m have also been reported [3]. The level of conservation is particularly striking given the high mutation rates seen in RNA viruses.
The s2m sequence is also present in the genome of the severe acute respiratory syndrome (SARS) coronavirus (SARS-CoV) where its three-dimensional crystal structure has been characterized in great detail [4]. The function of s2m remains obscure, although in the case of SARS-CoV, it was suggested that the structure might interfere with protein synthesis through mimicry of small subunit ribosomal RNA (SSU rRNA) and subsequent binding of ribosomal proteins [4]. This proposed affinity was, however, only observed when interactions were modeled using prokaryotic proteins. One of the many interesting features of s2m is that there seems to be a high level of conservation on all levels; primary structure (sequence) of both stem and loop regions, secondary (stem/loop structure) and tertiary conformation, indicating that all of these characteristics are important for functionality. The conserved nature of s2m has led researchers to suggest using it as a drug target [4] and there are also protocols describing how primers targeting the conserved s2m sequence can be used for virus discovery using reverse transcription PCR [5].
The phylogenetic distribution of s2m seems to support a model where the genetic element can be transferred horizontally and mobility of s2m has been suggested in the literature [3]. We present a thorough analysis of the distribution of s2m in viral genomes and perform likelihood-based phylogenetic analyses of all the relevant virus groups. We propose models for the evolutionary history of s2m within the different host virus families and perform likelihood ratio tests to investigate the apparent mobile nature of this sequence motif.

Results
A total of 682 s2m-containing sequences were found in GenBank when a consensus sequence-based approach was used directly. Allowing a single nucleotide mismatch increased this number to 702 and when two mismatches were allowed, a total of 706 sequences could be found from four different virus families (Additional file 1). Within this set of accessions, representatives from all the s2m-containing virus families were found that contained two copies of s2m (Table 1), but no sequences were found to contain three or more. The s2m sequences from genomes with two copies were never identical to each other. The genomes of dog norovirus strains GVI.1/HKU_Ca026F/2007/HKG and GVI.1/ HKU_Ca035F/2007/HKG had identical 3'ends, both with two copies of s2m, although there were minor sequence differences in the rest of the sequences.
For the astroviruses and caliciviruses, partial RNAdependent RNA polymerase (RdRP) amino acid sequences were used for the phylogenetic analyses (257 and 265 amino acid residues, respectively), whereas larger parts of the polyprotein sequence could be unambiguously aligned for the picornaviruses (998 residues) and the coronaviruses (2833 residues). The s2mcontaining viruses did not form monophyletic groups in any of the four virus families (Figures 1A, 2, 3 and 4).
Except for the coronaviruses, numerous sequences were obtained from isolates that had not been assigned a specific taxonomic placement within their respective families in the NCBI Taxonomy database, but in general, all trees constructed were consistent with previously published reports on viral phylogeny.
The ln likelihood difference (Δln) between unconstrained optimal trees and the optimal tree where all s2m sequences were required to be holophyletic (i.e. a single origin of the s2m sequence) ranged from coronaviruses (3102.9), picornaviruses (939.0), astroviruses (392.3), to caliciviruses (3.4). In all but the case of the caliciviruses, monophyly of the s2m-containing viruses was rejected with p < 0.01 (Additional file 2: Figure S1). Furthermore, even constrained trees assuming two gains of s2m were strongly rejected for coronaviruses, picornaviruses and astroviruses (Additional file 2: Figure S2, S3 and S4). For example in the picornaviruses one explanation for the distribution of the s2m sequence could be that the sequence was gained four times ( Figure 3). The Approximately Unbiased (AU) test found that the tree where the s2m sequences were constrained into a single clade (requiring holophyly of s2m containing taxa) was significantly worse than the optimal tree (Additional file 2: Figure S3). Even constraining Bat picornavirus 3 and Canine picornavirus into a single clade was significantly worse than the optimal tree (Additional file 2: Figure S3).

Discussion
Traditionally, the term mobile genetic elements (MGEs) has been restricted to include bacteriophages, plasmids and transposons, although it is now widely recognized that this classification is becoming obsolete as many elements with novel features as well as new combinations of known features are found [7,8]. Homologous and non-homologous recombination events have been inferred based on observations done in a wide range of  virus groups, particularly in the single-stranded RNA viruses [9]. The genomes of bacteriophages are also known to contain genetic material from multiple sources and appear to be quite promiscuous when it comes to acquiring novel genomic features [10]. For non-phage viruses, mobile genetic elements seem very rare. To our knowledge, the only other example of what appears to be a mobile genetic element in regular viruses is the S7 domain found in certain double-stranded RNA viruses, although the degree of amino acid conservation is quite low [11].
No highly supported branches separated the dog norovirus sequences in the calicivirus group ( Figure 1A), so a parsimonious model for the distribution of s2m here would be a single gain of this character. In contrast to what was found when investigating the other three s2mcontaining virus families, a single gain of s2m was not rejected by the AU test. Phylogenetic analyses of the s2m sequences themselves ( Figure 1B) could not, however, resolve the evolutionary history of s2m for this group, most likely due to residues either being too rapidly evolving or too conserved to give a good phylogenetic signal (data not shown). For the coronaviruses, two gains and two losses would have to be postulated in order to explain the distribution of s2m in the most parsimonious way (Figure 2). If the ancestral state for the gamma/delta coronavirus group was to contain s2m, losses would have to be proposed for the Night heron and Wigeon isolates, albeit there was little support separating these species so they could represent a single loss. Loss would have to be proposed, however, for the Beluga whale coronavirus. The s2mcontaining, monophyletic group comprising the SARS virus and a bat coronavirus can be explained through a single gain. In the picornaviruses, there appears to a more complex distribution of s2m ( Figure 3). Single gains can explain the two monophyletic groups that have s2m (paraturdiviruses and Equine rhinitis virus B 1 and 2) and the presence of s2m in Pigeon picornavirus B, but the phylogenetic placement of the two other s2m viruses (Bat picornavirus 3 and Canine picornavirus) is more ambiguous. They are not separated by highly supported branches and could thus reflect a single gain, or more complex explanations can be proposed, implying multiple gains and losses.
The broad distribution of s2m in the astroviruses (Figure 4) has been noticed previously [3] and a possible explanation could be that the ancestral state for this entire family was to contain s2m. For this to be true, a single loss would have to be proposed for the avastroviruses (Avastrovirus 3 and Astrovirus strain CDB-2012). Two losses would have to have occurred in Bat astrovirus 1 and Bat astrovirus strain Tm/Guangxi/ LD38/2007. All members of the large, monophyletic group that includes the classical human astroviruses and has the California sea lion astrovirus 2 as the most basal branch contain s2m. This group is not separated from the rest of the (non-s2m containing) mamastroviruses by any highly supported branches, and is thus possible that this is a 'primitive' member of this virus family and that the absence of s2m in the remaining isolates can be explained through a single loss.
It is intuitive that loss of a complex character, and in particular a character that can only provide an evolutionary advantage in a direct or indirect interplay with an existing cellular mechanism and does not lead to a tremendous increase in fitness, is more likely than gain of such a feature in an evolutionary perspective. Given the data currently available for s2m we conclude that it is impossible to establish a statistical model that can take into account any such differences. In spite of this, we still believe that horizontal transfer is the most plausible explanation for the distribution of s2m. It is also formally possible that this is a case of convergent evolution, but we consider this highly unlikely given the high degree of s2m similarity and the complexity of the character. Alternative hypotheses would have to propose that s2m was present in the last common ancestor of the + ssRNA viruses. There is also an apparent lack of intermediate/primitive forms of s2m motifs in GenBank. It is probable that our search strategy for s2m(−like) motifs has a certain false negative rate, but the number of sequences that could be found quickly reached a plateau, where allowing more substitutions or fixing a smaller number of consensus motif characters only led to an exponential increase in the number of obvious false positives (data not shown).
A model where containing the s2m motif is the proposed ancestral state for all + ssRNA viruses would have to postulate a large number of independent losses and a extraordinary selection pressure to maintain s2m in certain viral lineages. This seems unlikely as s2m only appears to provide a somewhat subtle (yet immediate) selective advantage for the host viruses. The fact that s2m remains conserved in spite of the high mutation rates seen in RNA viruses indicates that the virus somehow benefits from acquiring the sequence motif, but unfortunately there are very few examples from the published literature on closely related viruses that differ in their s2m status. The turdiviruses (Figure 3) were all collected from dead birds, and the authors were unsuccessful in their efforts to culture the viral strains for further characterization [12]. All bats that were found to contain bat picornaviruses appeared healthy (Figure 3 and [13]) and there is no data indicating that the two viral strains that apparently have lost s2m within the delta/gammacoronavirus group (Night-heron coronavirus strain HKU19-6918 and Wigeon coronavirus strain HKU20-9243) are significantly different from the other members of this group in terms of pathogenicity, host specificity, etc. [14]. Exchanging the (non-s2m containing) 3'-end of a murine coronavirus (MCV) with the (s2m containing) 3'-end of a SARS-CoV did not appear to have a dramatic effect on the virus [15] and an IBV strain with a deleted version of s2m that was discovered as an escape mutant in a vaccine development project (GenBank accession number JF274479) did not appear phenotypically different from closely related viruses in culturing experiments (Dr. Shengwang Liu, personal communication).
The viruses that contain s2m can infect a wide range of higher vertebrates, including birds, bats, horses, dogs and humans, and display different tissue tropisms. The most likely scenario for the emergence of a new s2mcontaining virus would be a situation where a co-infection includes both an s2m-containing donor virus and a recipient virus. Bats have been shown to carry many different viruses, including members of all the families that have been shown to harbor s2m [13,16]. Due to their mobility, feeding habits, long life span, roosting behavior, general virus susceptibility, etc., it has been proposed that bats may represent an important reservoir for emerging viruses [17]. However, for the coronaviruses, several s2mcontaing members have been postulated to have an avian origin [14]. On a molecular level, is has been suggested that transfer of s2m occurs through non-homologous recombination in an replication-dependent manner [3]. Based on our sequence alignments, it is impossible to determine whether or not such a model should include just the hairpin structure or if transfer of the entire 3' end of the genome represents a more likely scenario. The noncoding nature of this part of the genome is associated with high mutation rates and any sequence similarity in the s2m-flanking regions (particularly downstream, near the poly(A) tail) would quickly be lost due to the high error rates observed in replicating RNA virus genomes [18].
In most cases where phylogenetic analyses are used to investigate the horizontal transfer of a genetic element or a DNA-containing organelle such as plastids or mitochondria, it is possible to do a phylogenetic analysis of the genetic element that has (presumably) been transferred and then compare the resulting topology with that of the hosts. s2m is short, but due to its secondary structure it can be unambiguously aligned and there should also be sufficient characters that show some degree of variability to give reasonable resolution if data from closely related species are compared. Regardless of this, we were unable to find any correlation between s2m mutational patterns and host phylogenies. For instance, we were unable to assess whether the two copies of s2m found in the dog norovirus strains came from independent sources or if they are the result of some sort of duplication event and subsequent independent evolution ( Figure 1B and Table 1), and the phylogeny of the deltacoronovirus s2m sequences was poorly resolved and did not match the phylogeny of the hosts (data not shown). There seems to be a number of loci that 'permit' certain substitutions and that these mutate quite quickly, masking any phylogenetic signal.
Although our analyses did not address the function of s2m per se, we believe that our observations might provide some clues as to how s2m might evolve and provide a selective advantage to the host viruses. We believe that s2m must have some sort of 'autonomous' function that does not require complex interactions with other parts of the viral genome/transcriptome as there do not appear to be any conserved flanking regions in s2m viruses (proximity to coding region, adjacent reading frames, nucleotide motifs etc.). Neither do there appear to be any conserved amino acid motifs in any of the annotated open reading frames except the GDD core of the RNAdependent RNA polymerase when looking at protein sequence data from representatives from the four virus families (data not shown). The conserved nature of s2m might also imply that the s2m targets are homologous as all infected organisms are (relatively) closely related in an evolutionary perspective whereas the viruses are distantly related when looking both at molecular data and functionality (genome replication and transcription/ translation strategies). An intracellular target also seems plausible, as all s2m containing viruses replicate in the cytoplasm and this is where s2m is likely to be available for interactions with the cellular machinery. The observation that s2m can apparently be transferred between unrelated viruses and remain functional (under selection pressure to maintain sequence and structure) also suggests strongly that the target for s2m is host-specific and not viral. Based on structural similarities between micro RNA (miRNA) hairpins involved in gene regulation, we propose that s2m functions through a RNA interference (RNAi)-like mechanism, possibly targeting homologous sequence loci in infected organisms. Recent observations using a reverse genetics-based approach and a recombinant Sindbis virus indicate that the required cellular machinery for this to function should be in place in human cells [19], and this model would also be consistent with an additive effect, where more copies of s2m would allow the formation of more miRNA/protein complexes, resulting in a more profound effect on target gene regulation.
Based on the finding of s2m in what appears to be newly emerging viruses, such as SARS-CoV, we also believe that s2m still maintains its mobility and will play a role in the future of virus evolution.

Conclusions
The s2m sequence motif appears to be an active mobile genetic element that thus far can be found in four different families of + ssRNA viruses. It seems likely that s2m provides some kind of selective advantage for the viruses that contain the motif, and a possible function could be related to RNAi-like gene regulation of infected organisms.

Methods
All sequence data were downloaded from GenBank. Only accessions where continuous sequence information from the 3' end of the genome to the locus selected for phylogenetic analyses was available were used when investigating the evolutionary relationship between the (s2m-containing) viruses.
By using all the available sequence information from publications pertaining to s2m, several conserved sequences domains within s2m could be identified. These sequence motifs were used to perform nucleotide BLAST searches through the NCBI portal and >400 s2mcontaining sequences could be identified. The s2m motifs (43-44 nucleotides long) were individually extracted and aligned, and the following consensus sequence could be generated: CGNGG(N)CCACGNNGNGT(N)ANNANCG AGGGT(N)ACAG (N's in parenthesis indicate possible indels) for the conserved core region of s2m. This text string profile was used to search all viral sequences in GenBank using different combinations of the indicated indels and allowing for nucleotide substitutions.
All the alignments were constructed using the Clustal W algorithm [20] and manually edited using Bioedit (version 7.0.5.3; [21]). Only unambiguously aligned domains were included in the subsequent analyses. Phylogenetic analyses were performed using RAxML with the optimal GTR (General Time Reversible) with Gamma distribution model for amino acid substitution and 100 bootstrap replicates [22]. Minimally constrained trees were constructed to test whether trees with monophyly of s2m-containing sequences were significantly less likely than the optimal unconstrained tree. The primary constraint was to make all s2m-containing taxa monophyletic. If more than two s2m-containing clades were present alternate topologies were constructed in a pairwise manner to test topologies with two s2m clades constrained together. The optimal tree compatible with the constraint was calculated using RAxML with the above amino acid substitution model, followed by site likelihood calculation with RAxML. The Approximately Unbiased (AU) test was calculated using CONSEL [23].

Additional file
Additional file 1: List of s2m-containing sequences found in Genbank.
Additional file 2: Figure S1, S2, S3 and S4. For each virus family, the optimal tree with bootstrap support values is shown (identical to Figures 1, 2, 3 and 4). Each of the s2m-containing clades was numbered, from top to bottom without regard to taxonomy, and the total number of observed s2m gains was mapped onto the tree. Constraint trees were constructed to test a series of hypotheses for s2m gains. For example, all s2m-containing sequences were constrained to a single clade to see if a tree with only a single gain of the s2m sequence was significantly less likely given the alignment. Other less constrained trees were also constructed to test if trees with more than one s2m gain were also significantly less likely given the alignment. The specific constraints on s2m clades are shown in the table below. In all cases the differences in log likelihood and p-values for the different trees are shown in comparison to the most likely unconstrained tree. In addition to the Approximately Unbiased (AU) p-values, other test support values are also shown (np = bootstrap probablility, bp = bootstrap proportion, pp = posterior probability, kh = Kishino-Hasegawa, sh = Shimodaira-Hasegawa, wkh = weighted Kishino-Hasegawa, wsh = weighted Shimodaira-Hasegawa).