dataesr/scientific-paragraphs-categorization huggingface dataset:
The RMIP values of 0. 05 or greater were considered signi cant  44. Genomic regions identi ed by both methods-meaning low empirical p-values and high RMIPs-are high-con dence hits; those identi ed by only one are still signi cant but of lower con dence. For those traits that did not pass through signi cant thresholds and also showed weak signals in the above two methods, GWAS was also performed with FarmCPU in GAPIT package from R software without any iterations. GWAS analysis identi ed a total of 200 MTAs for different aphid resistance-related traits, owering time, and plant height which are discussed below (Supplementary Table  S2 -a). AC1-Aphid Count rst rating dates and AC2-Aphid Count second rating dates FarmCPU-RMIP model selection analysis identi ed 11 loci for AC1(Fig. 1 , Supplementary Table  S2-a ).
Of these, eight genes encoding leucine-rich repeats (LRR) were found on SBI-05, SBI-06, and SBI-07. The SNP marker, S6_334458, associated with aphid damage (D1), is 5 kb away from Sobic. 006G001900 encoding a LRR gene that was upregulated during aphid infestation in the previous study  45. Four genes near NDRE associated SNPs on SBI-03, SBI-06 and SBI-08 were similar to lipoxygenase (LOXs), and two genes on SBI-05 and SBI-07 encoded avonol biosynthesis genes (Supplementary Table  S4 ). The NDVI marker, S1-7442679, was found within ~99 kb of the Sobic. 001G0955 gene that encodes a WRKY1 transcription factor. Additionally, we pooled gene datasets obtained from Tetreault et al. 46  that were detected within a 200 kb distance of MTAs. We analyzed genes that were induced in both resistant and susceptible plants upon sugarcane aphid infestation using 16 gene-expression modules in this study.
46  (Supplementary Table  S4  and S5 ). When all the modules from the Tetreault et al. 46  study were used for candidate gene analysis, one 12-oxophytodienoic acid reductase gene (Sobic. 004G180500) was found very close within 23 kb from S4_53265404. The Sobic. 007G151100 encoding 12oxophytodienoic acid reductase was also found within 483 kb to a marker, S7_58778565, associated with aphid count (AC1). We also found three LRR genes on SBI-05 and SBI-07, and one WRKY1 gene on SBI-01, which were commonly found in the studies by Kiani and Szczepaniec, 45 and Tetreault et al. 46. We also analyzed the number of genes and gene frequency on individual sorghum chromosomes found near our SNPs associated with different aphid traits, by using these two datasets  45, 46. The highest gene frequency distribution was observed on SBI-01 followed by SBI-06 (Supplementary Fig.
When we compared our markers with these FSt outlier markers for sugarcane aphid resistance, 34 markers were found to be in close proximity (within 500 kb) to 40 markers identi ed in Muleta et al. 56  (Supplementary Table  S7 ). These markers were as close as 0 to 356 kb between the two studies and most of these markers came from SBI-06-and SBI-09. There are only a few studies that have elucidated the changes in gene expression upon sugarcane aphid infestation in sorghum. In our study, we were able to detect 6% of known unique genes from these two datasets related to sugarcane aphid response using our MTAs identi ed. Our candidate gene search was limited to these two datasets from the existing two experiments, and our window to examine these genes was only 200 kb within signi cant markers  45, 46.
High throughput phenotyping using UAS The remote sensing datasets used in this study consisted of multispectral imagery data collected over sorghum eld trials in Tifton, GA, during six eld visits (six ight dates) from July through September 2020. Only the data from Aug. 25 were used in this study as it coincided with visual scoring dates and had higher correlation (Supplementary Table  S1 -d). The imagery was collected with a Micasense RedEdge 5band sensor using a DJI Matrice M100 UAS. Spectral bands collected were blue (B, 475 +/-20nm), green (G, 560 +/-20nm), red (R, 668 +/-10nm), red edge (RE, 717 +/-10nm), and near infrared (NIR, 842 +/-40nm). To georectify the image mosaics, ground control points (GCPs) were established at the beginning of the study and marked with a temporary landmark, using a Trimble Geo7X global navigation satellite system receiver with real-time kinematic correction (i.
e. , RTK), providing +/-2 cm horizontal accuracy at the coordinate location. Prior to each ight, the GCPs were targeted, and calibration images were collected using Micasense calibration protocols. Pix4DMapper software (v. 4. 6. 4; Pix4D SA, Switzerland) was used to produce orthorecti ed mosaics, photogrammetric point clouds and surface elevation models, georeferenced to within a satisfactory root mean square error (RMSE). The nal spatial resolution of the mosaicked images was 2. 64 cm. The UAS was own with an iPad Mini using Pix4DCapture ight control software, and ight parameters were ~38 meters altitude, at 4. 92 m s -1 (11 mi h -1 ), resulting in 80% side and front overlap in the RedEdge camera imagery. The three vegetation indices calculated for this study were Normalized Difference Vegetation Index (NDVI), Normalized Difference Red-Edge (NDRE), and Soil Adjusted Vegetation Index (SAVI).
GWAS with FarmCPU-Model selection RMIP-method To complement the GLM analysis above, model-selection GWAS was performed with FarmCPU 42 as part of the GAPIT3  91, 92  software package for R. The same genotype le was used for all traits, and the model included three genetic principal coordinates and a kinship matrix calculated using the VanRaden method  93  with the function GAPIT. kinship. VanRaden(). Early tests found that the FarmCPU results, like many model selection algorithms  94  , were unstable when six samples were accidentally excluded from the analysis and generated very different lists of signi cant loci. To overcome this instability, we performed model selection in two steps. First, a null distribution of results created by randomly permuting the genotypes relative to the phenotypes 100 times.
Phenotypic Variance Explained To estimate variance explained by each marker (or set of markers), we used GAPIT3  91, 92  to calculate a kinship matrix using the VanRaden method  93  and the prcomp() function in R to derive genetic principal components (PCs). For each phenotype, we t a model including the rst three PCs, the kinship matrix, and any covariate phenotypes (e. g. , owering time), as with the GWAS analyses. The residuals from that model were then t in a second linear regression with the SNP(s) of interest to determine the fraction of residual variance explained in the form of the model R 2. Code for this analysis is available on GitHub at https://github. com/wallacelab/paper-punnuri-aphids-2021.
GWAS with single-run FarmCPU For those traits that did not pass-through signi cant thresholds and also showed weak signals in the above two methods, GWAS was also performed with FarmCPU in GAPIT package from R software without any iterations. The results from this analysis are displayed in the supplementary les and were only performed in the case of aphid count second-dates rating (AC2) and greenhouse evaluation data as they showed genes relevant to sugarcane aphid herbivory on sorghum. The results need to be further validated as no iterations were conducted in this method. Candidate gene search from previous gene studies Candidate genes were selected based on data available on differentially expressed genes (DEGs) between resistant and susceptible sorghum from two previous experiments  45,.
Acknowledgements The authors would like to thank  Trevor Perla ,  Michael Purvis , and  Tyler Bailey  for technical support, as well as  Barb Driskel ,  Hannah Baker  and  Cheyenne Edwards  for running the greenhouse phenotyping trials. We would like to thank the  USDA NIFA Capacity Building grant program , award no  2019-38821-29057  for funding support. We would like to thank  USDA, Agricultural Research Service , National Plant Germplasm System and GRIN resources for providing SAP seeds continuously throughout the study period. Also, we would like to thank  Dr. Xinzhi Ni  for providing seeds for the rst year of the SAP planting. The research was also supported in part by the  Non-Assistance Cooperative Agreement   58-6048-6-016 , which was funded by the  United States Department of Agriculture -Agricultural Research Service  with the research project entitled " Areawide Pest Management of Invasive Sugarcane Aphid in Grain Sorghum  ( 62048-21220-018-00D ).
Results and discussion Co-phylogenetic history of an insect-plant antagonistic interaction. First, we created an extensive phylogenetic dataset including seven genetic markers for 71% of swallowtail species diversity (408 of ~570 described species, see 'Methods'). This dataset leads to the assembly of the most complete and wellresolved dated phylogeny of swallowtail butterflies (79% of nodes with strong bootstrap support defined as ≥95%; Supplementary Figs. 1 2 3 ). Both tribe-and genus-level relationships are mostly consistent with previous results using multilocus datasets  18, 34, 35, [42] [43] [44] [45]. However, our species tree benefits from a phylogenomic backbone that we recently inferred at the genus level for the Papilionidae using genome-scale data  46. Second, we compiled host-plant preferences for each swallowtail species in the dataset, and we performed ancestral-state estimations (see an ancient and highly conserved association with Aristolochiaceae throughout swallowtail butterflies evolution.
Reconstructions of co-phylogenetic history for other insect-plant antagonistic interactions have shown either synchronous diversification  11  or herbivore diversification lagging behind that of their host plants  10, 47. We assembled a molecular dataset for ~49% of the species diversity of Aristolochiaceae (247 of ~502 described species; see 'Methods') and reconstructed their phylogeny (Supplementary Fig. 7 ), which is in agreement with previous works  [48] [49] [50] [51] [52]. Divergence time estimates strongly suggest synchronous radiations of Papilionidae (55. 4 million years ago [Ma], 95% credibility intervals (CIs): 47. 8-71. 0 Ma) and Aristolochia (55. 5 Ma, 95% CIs: 39. 2-72. 8 Ma) since the early Eocene (Fig. 2  and Supplementary Figs. 3, 8 and 9). This result is robust to known biases in inferring divergence times, with slightly older ages Fig.
To characterize the macroevolutionary patterns of the Aristolochia/Papilionidae arms race in space, we assembled two datasets of current geographic distributions for all species included in the phylogenies of both Aristolochiaceae and Papilionidae. We reconstructed the historical biogeography of both groups, taking into account palaeogeographical events throughout the Cenozoic (see 'Methods'). Along with the known fossil record of both groups  [53] [54] [55] [56] [57]  , these results suggest that both Papilionidae and Aristolochia were ancestrally co-distributed throughout a region, including West Nearctic (WN), East Palaearctic (EP), and Central America (CA) in the early Eocene, when Asia and North America were connected by the Bering land bridge (Fig. 2  and Supplementary Figs. 10  and 11 ).
This supports the hypothesis that insect-plant interactions contributed to P450-gene family diversification, with P450s being key innovations that explain the evolutionary and ecological success of phytophagous insects  14, 15, 36, 38, 68, 69. However, host-plant shifts not only alter single genes but may also influence unlinked genes  40. Moreover, host-plant shifts can accompany changes of the abiotic environment, which may, in turn, require further biotic adaptation (new predators and/or competitors). But the macroevolutionary and genomic consequences of the evolutionary dynamics of host-plant shifts have not yet been demonstrated. Relying on a genomic dataset comprising 45 genomes covering all swallowtail genera  46, [70] [71] [72]  , we constructed two specific datasets (Dataset 1: 520 genes and Dataset 2: 1533 genes; mean gene coverage = 26.
Using a maximum-likelihood (ML) method, we estimated the ratio of non-synonymous substitutions (dN) other synonymous substitutions (dS) in all branches where a host-plant shift was identified relative to branches with no host-plant shift  75, 76  (see 'Methods'). The dN/dS analyses on branches with host-plant shifts (combined or not with environmental shifts) showed more genes with a subset of codons evolving under positive selection (dN/dS > 1) in lineages shifting to a new plant family, although the difference was marginally non-significant for the smallest dataset and highly significant for the second dataset containing more genes (Fig. 4b , Supplementary Fig. S19  and Supplementary Fig. 3  Host-plant shifts lead to repeated bursts in diversification rates and a sustained overall increase in diversification through time.
2 , P = 0. 0501/0. 0079 for the two datasets, respectively, Wilcoxon rank-sum test, see 'Methods' for the definition of the datasets). However, dN/dS analyses on branches with environmental shifts indicated a balanced number of genes under positive selection (Fig. 4c , Supplementary Fig. S19  and Supplementary Table  2 , P = 0. 336/0. 8162 for the two datasets, respectively, Wilcoxon rank-sum test), suggesting a lower impact of environmental shifts than host-plant shifts. We then performed dN/dS analyses for branches with host-plant shifts only (not followed by environmental shifts) and found that swallowtail lineages shifting to a new host-plant family had significantly more genes under positive selection (4. 41%/3. 98% of genes under positive selection for the two datasets, respectively; Supplementary Table  2 ) than non-shifting lineages (3.
02%/2. 43% of genes under positive selection for the two datasets, respectively, Fig. 4d , Supplementary Fig. S19  and Supplementary Table  2 , P = 0. 0071/0. 00156 for the two datasets, respectively, Wilcoxon rank-sum test). Surprisingly, the dual changes in environment and host-plant preferences did not spur molecular adaptation across swallowtail lineages compared to control branches (P = 1/0. 4439 for the two datasets, respectively, Wilcoxon rank-sum test; Fig. 4d , Supplementary Fig. S19  and Supplementary Table  2 ). Comparing the proportion of genes under positive selection between the branches with dual changes and branches with host-plant shifts only shows a marginally significant difference with Dataset 1 and no difference with Dataset 2 (P = 0. 0327/0. 1471 for the two datasets, respectively, Wilcoxon rank-sum test; Fig.
4d  and Supplementary Fig. S19 ). However, this result might be an artefact due to the use of a few branches to perform the statistical comparison. Although we did not control for the effect of multi-nucleotide mutations  77  , which should affect dN/dS analyses equally for control and host-plant shift branches, we checked individually the gene alignments and performed sensitivity analyses, which showed that our results are not driven either by an excess of misaligned regions or missing data and GC-content variations among species (see 'Methods' and Supplementary Figs. 20-26). Finally, given that fixing the topology for CodeML (see 'Methods') can spuriously inflate substitution rates on some branches  78  , we computed the proportion of genes under positive selection by selecting the gene trees from the largest dataset (Dataset 2) for which the focal branches were recovered (in agreement with the species tree).
These analyses confirmed the previous results suggesting more genes under positive selection during host-plant shifts (P = 0. 0444, Wilcoxon rank test; Supplementary Table  2 ). We further studied the functional categories of positively selected genes by using gene ontology (GO) analyses (PANTHER and EggNOG; see 'Methods'). Applied to the high-quality genomes of Papilio xuthus  71  and Heliconius melpomene  79  , we found that ~70% of the genes could be associated with a gene function and ~30% lacked annotation, which suggests a gap of knowledge in the current insect database of gene function. Among the annotated genes, we found that genes under positive selection along branches with host shifts did not contain over-or under-represented functional GO categories: 252 out of 1213 GO categories represented by genes under positive selection (P > 0.
Methods Time-calibrated phylogeny of Papilionidae. We assembled a supermatrix dataset with available data extracted from GenBank as of May 2017 (most of which has been generated by our research group), using five mitochondrial genes (COI, COII, ND1, ND5 and rRNA 16S) and two nuclear markers (EF-1a and Wg) for 408 Papilionidae species (~71% of the total species diversity) and 20 outgroup species. We aligned the DNA sequences for each gene using MAFFT 7. 110  84  with default settings (E-INS-i algorithm), and the alignments were checked for codon stops and eventually refined by eye with Mesquite 3. 1 (available at: www. mesquiteproject. org). The best-fit partitioning schemes and substitution models for phylogenetic analyses were determined with PartitionFinder 2. 1. 1  85  using the greedy search algorithm and the Bayesian Information Criterion.
Following the recommendation of IQ-TREE developers, we also set smaller perturbation strength (-pers 0. 2) and a larger number of stop iterations (-nstop 500) to avoid local optima. We performed 2000 ultrafast bootstrap replicates to investigate nodal support across the topology, considering values ≥95 as strongly supported nodes  89. Estimating phylogenetic relationships for such a dataset is computationally intensive with Bayesian inference. The ML tree inferred with IQ-TREE was used as a starting tree for Bayesian inference as implemented in MrBayes 3. 2. 6  90. Rather than using a single substitution model per molecular partition, we sampled across the entire substitution-model space  91  using reversible-jump Markov Chain Monte Carlo (rj-MCMC). Two independent analyses with one cold chain and seven heated chains, each run for 50 million generations, sampled every 5000 generations.
After enforcing the fossil calibrations, we set the following settings and priors: a partitioned dataset (after the best-fitting PartitionFinder scheme) was analysed using the uncorrelated log-normal distribution clock model, with the mean set to a uniform prior between 0 and 1, and an exponential prior (lambda = 0. 333) for the standard deviation. The branching process prior was set to a birth-death 107 process, using the following uniform priors: the birth-death mean growth rate ranged between 0 and 10 with a starting value at 0. 1, and the birth-death relative death rate ranged between 0 and 1 (starting value = 0. 5). We performed four independent BEAST analyses for 100 million generations, sampled every 10,000th, resulting in 10,000 samples in the posterior distribution, of which the first 2500 samples were discarded as burn-in.
78-113. 70 Ma]  115  ). We used the estimated molecular ages of the different host-plant groups to constrain our inferences of ancestral host plants a posteriori. We preferred such an approach compared to a more constrained one in which the DEC model is informed with a matrix of host-plant appearances based on their estimated ages by implementing matrices of presence/absence of the character states through time (equivalent to the time-stratified palaeogeographic model, see below for inference of biogeographical history). Finally, the Bayesian approach implemented in BayesTraits 3. 0. 1  110  was performed to provide a cross-validation of ML analyses. This approach automatically detects shifts in rates of evolution for multistate data using rj-MCMC. The number of parameters and priors was set by default.
W. , July 2019) favour splitting them into four families: Aristolochiaceae (Aristolochia and Thottea), Asaraceae (Asarum and Saruma), Hydnoraceae (Hydnora and Prosopanche), and Lactoridaceae (Lactoris), collectively called the perianthbearing Piperales. Therefore, we extracted and assembled a supermatrix dataset with available data from GenBank for the perianth-bearing Piperales and its sister lineage, the perianth-less Piperales including Saururaceae and Piperaceae (as of May 2017, most of which has been generated by our research group). We obtained four chloroplast genes (matK, rbcl, trnL and trnL-trnF) and one nuclear marker (ITS) for 247 species of perianth-bearing Piperales (~49% of the total species diversity  130  ) and six outgroups from perianth-less Piperales. We could not include the two genera Hydnora and Prosopanche (Hydnoraceae) because available genetic data do not overlap those of perianth-bearing Piperales  126, 128, 131, 132.
We applied the same analytical procedure that we did for Papilionidae. DNA sequences for each gene were aligned using MAFFT 7. 110  84  with default settings (E-INS-i algorithm and Q-INS-I to take into account secondary structure). Resulting alignments were checked for codon stops and eventually refined by eye with Mesquite 3. 1 (available at: www. mesquiteproject. org). The best-fit partitioning schemes and substitution models for phylogenetic analyses were determined with PartitionFinder 2. 1. 1  85. All gene alignments were concatenated into a supermatrix; the final dataset is available in Figshare (see Data availability). Phylogenetic relationships were estimated with Bayesian inference as implemented in MrBayes 3. 2. 6  90. Rather than using a single substitution model per molecular partition, we sampled across the entire substitution-model space  91  using rj-MCMC.
We set the following settings and priors: a partitioned dataset (after the bestfitting PartitionFinder scheme) was analysed using the uncorrelated log-normal distribution clock model, with the mean set to a uniform prior between 0 and 1, and an exponential prior (lambda = 0. 333) for the standard deviation. The branching process prior was set to a birth-death 107 process, using the following uniform priors: the birth-death mean growth rate ranged between 0 and 10 with a starting value at 0. 1, and the birth-death relative death rate ranged between 0 and 1 (starting value = 0. 5). We performed four independent BEAST analyses for 100 million generations, sampled every 10,000th, resulting in 10,000 samples in the posterior distribution of which the first 2500 samples were discarded as burn-in.
All analyses were performed on the CIPRES Science Gateway computer cluster  94  , using BEAGLE  95. Convergence and performance of each MCMC run were evaluated using Tracer 1. 7. 1  92  and the ESS for each parameter. We combined the four runs using LogCombiner 1. 8. 4  97. The MCC tree was reconstructed with median age and 95% CI. The BEAST files generated for this study are available in Figshare (see Data availability). Dual biogeographic history of Papilionidae and Aristolochiaceae. We estimated the ancestral area of origin and geographic range evolution for both clades using the ML approach of DEC model  109  as implemented in the C++ version  137, 138  that is available at: https://github. com/champost/DECX. To infer the biogeographic history of a clade, DEC requires a time-calibrated tree, the current distribution of each species for a set of geographic areas, and a time-stratified geographic model that is represented by connectivity matrices for specified time intervals spanning the entire evolutionary history of the group.