s2fieldsofstudy: Biology Medicine

text:
Applying dynamic contrast enhanced MSOT imaging to intratumoral pharmacokinetic modeling

Examining the dynamics of an agent in the tumor microenvironment can offer critical insights to the influx rate and accumulation of the agent. Intratumoral kinetic characterization in the in vivo setting can further elicudate distribution patterns and tumor microenvironment. Dynamic contrast-enhanced Multispectral Optoacoustic Tomographic imaging (DCE-MSOT) acquires serial MSOT images with the administration of an exogenous contrast agent over time. We tracked the dynamics of a tumor-targeted contrast agent, HypoxiSense 680 (HS680), in breast xenograft mouse models using MSOT. Arterial input function (AIF) approach with MSOT imaging allowed for tracking HS680 dynamics within the mouse. The optoacoustic signal for HS680 was quantified using the ROI function in the ViewMSOT software. A two-compartment pharmacokinetics (PK) model constructed in MATLAB to fit rate parameters. The contrast influx (kin) and outflux (kout) rate constants predicted are kin = 1.96 × 10−2 s-1 and kout = 9.5 × 10-3 s-1 (R = 0.9945).

Introduction
Multispectral optoacoustic tomography (MSOT) is a novel in vivo imaging modality with high spatial and temporal resolution at clinically significant detection depths [1][2][3]. Near-infrared (NIR) light pulses are absorbed by chromophores in the tissue, leading to heat production and thermoelastic expansion. This results in acoustic signal production which can be detected by an ultrasound transducer. MSOT imaging takes advantage of the optoacoustic (OA) effect with ultrasound detection of optical excitation [4]. NIR light with wavelength between 600-900 nm do not attenuate in most tissues, resulting in high resolution images at improved depth than optical imaging [1]. The use of OA overcomes the main limitations of optical and ultrasound imaging including limited tissue penetration (< 8 mm), and lack of potential for molecular information using a wide array of contrast agents [5,6].
Because each chromophore absorbs light and emits a distinct optoacoustic signal, it is possible to detect multiple chromophores simultaneously in high-resolution images as well as to obtain biological information, such as oxygenation state and tissue inflammation [2,[7][8][9][10]. This ability to simultaneously track multiple chromophores based on absorption spectra allows for one multispectral scanning session to yield multiple images depicting the biodistribution of each chromophore [11][12][13]. Thus, endogenous and exogenous agents can be used in combination to facilitate the identification of tumors. One hallmark of tumors is a region of tissue hypoxia from abnormal microvessel formation and decreased diffusion, leading to lack of vascular flow [14]. HypoxiSense 680 (HS680) targets and binds to carbonic anhydrase IX, a marker for tumor hypoxia [15]. HS680 absorbs at 670 nm and emits at 685 nm ( Fig. 1), which can be tracked using MSOT imaging along with oxygenated hemoglobin (HbO 2 ).
In addition to detecting distinct spectral signatures, optoacoustic imaging is an intrinsically fast imaging technique [16]. High temporal resolution allows the advantage of visualizing dynamic, biological events, which is difficult in ultrasound imaging due to motion artifacts [17][18][19][20]. One application is a pharmacokinetics (PK) study. Pharmacokinetic modeling maps out the kinetics of an administered drug or contrast agent in the plasma of an organism, thus playing a predictable role in the dose-response relationship in toxicity assessment [21][22][23].
Here we used a 2-compartmental PK model to simplify the body of an organism into one or two compartments, a commonly used approach [24,25]. in vivo data are collected and fitted against a PK model with estimated rate parameters.
Combining the functionality of contrast agents with a fast tracking imaging modality yields a powerful tool that can elucidate contrast kinetics and agent biodistribution in vivo [18]. DCE-MSOT has not been used previously to investigate PK with a contrast agent targeting tumor hypoxia [18,26]. This study utilizes HS680-MSOT contrast-imaging to study intratumoral pharmacokinetics in the in vivo setting.

Materials
The powder form of HS680 (Perkin-Elmer Inc., Waltham, MA, USA) was diluted to a concentration of 6.1 mM in saline solution. The absorbance and fluorescent emission of HS680 were measured using a FP-8300 spectrofluorometer (Jasco, Oklahoma City, OK, USA) and the OA spectra was measured using an MSOT inVision 256-TF (iThera Medical GmbH, Munich, Germany) within a tissue mimicking phantom as in previous studies ( Fig. 1) [8,27]. A tail vein catheter was set up to inject each mouse with 100 μL of HS680 and a concentration of 4.5 μmol/kg mouse and followed by a catheter flush with 100 μL of saline solution. Each injection took ∼30 s.

Mouse breast tumor model
Three female, athymic mice, weighing 20 g each, were injected with 1 × 10 6 MDA-MB-231 cells in the mammary fat pad. After injection, mice were housed under pathogen-free conditions following the guidelines of the American Association for Laboratory Animal Care. When tumors were 7 mm in diameter (45 days post injection), PK and MSOT imaging studies were performed in accordance with the current standards of the University of Louisville Animal Care Facility, the Institutional Animal Care and Use Committee of the University of Louisville, the National Institutes of Health, the U.S. Department of Agriculture, and the U.S. Department of Health and Human Services.

Multispectral optoacoustic tomography (MSOT) imaging
A commercially available MSOT inVision 256-TF (iThera Medical GmbH, Munich, Germany) was used for optoacoustic measurements. The MSOT emitted near-infrared laser pulses at 10 Hz and wavelengths between 680 to 900 nm A 256 element, toroidal, ring-shaped transducer received optoacoustic signals. This 270°transducer setup allowed for cross-sectional tomographic imaging with an in-plane resolution of 75 μm which has been further confirmed in our previous manuscript [28], and a center frequency of 5 MHz. This arrangement provides clear visualization of the spatial distribution of injected agents in the animal, avoiding the "masking" of deeper pixels by contrast uptake in a shallow region. In addition, implantation of the orthotopic breast tumor in the mammary fat pad avoided high signal noise level from internal organs.
A mouse was anesthetized with 1.5% isoflurane and 0.8 l medical air with 0.1 l O 2 gas in preparation for imaging. A 27 gauge catheter was inserted in preparation for tail vein injection. To limit physiological changes in the animal and manual introduction of vessel volume fluctuation, we ensured that mouse physiology was stable for 10 min while the mouse was located within the MSOT unit before data collection and performed bolus injection over 30 s followed by a slow saline flush [29]. MSOT signal data were collected 5 min prior to injection of HS680, providing a baseline reading in the tumor region that represents the hemoglobin level. Multispectral images were acquired with three contiguous image slices centered on the tumor for 20 min following injection of HS680 using three wavelengths 680, 760, and 900 nm and 5 pulses per wavelength.
Advanced Molecular Imager (AMI-1000X) was used to detect the biodistribution of HS680 in the mouse model 24 h after injection.

Image reconstruction and spectral unmixing
A standard back projection algorithm from the ViewMSOT 3.6 software suite (iThera Medical) was used to reconstruct the images [30]. The linear spectral unmixing method was used to detect and separate unique signals for contrast agents (i.e. hemoglobin and HS680) [31,32] by comparing the optoacoustic signals measured at each pixel with the known absorption spectra of the target contrast agent(s), allowing for separation and visualization of the distribution of each agent within the animal. In the current study, this method was utilized to distinguish between signal differences in HS680 and hemoglobin in the tumor and blood vessels.

Pharmacokinetic modeling
The PK model for this study is a 2-compartmental PK algorithm custom-designed using MATLAB. The two compartments represent the plasma and the breast tumor microenvironment (Fig. 2). The mathematical approach is similar to the 2-compartmental models previously utilized in previous DCE-MRI and DCE-PET studies [33,34]. A differential equation is used to describe the kinetics of HS680 within the tumor microenvironment (Eq. (1)). At a given time point, C tumor (t) is the concentration (nM) of HS680 in the tumor microenvironment; C plasma (t) is the concentration of contrast in the plasma. The rate constant k in (s −1 ) describes the rate at which HS680 leaves the plasma and is taken up into the tumor space and k out describes the rate at which the agent returns from tumor space into the plasma. The differential equation assumes a known working concentration, which is 0.91 mM of HS680.
This in vivo experiment measures optoacoustic signal of HS680. The absorbance and fluorescent emission of HS680 was determined using a spectrofluorometer. (B) The optoacoustic spectrum was determined by inserting HS680 into a tissue mimicking phantom and measured using MSOT. The optoacoustic spectrum was utilized to identify HS680 in vivo as detected using MSOT.
MSOT imaging provides continuous OA signal acquisition at selected near-infrared wavelengths using HS680, acquiring OA signal for background and contrast shown in Eq. (3). Then, OA signal for HS680 is isolated and substituted with the HS680 concentration expression in Eq. (2). Isolation of C HS680 derives an expression for HS680 concentration in relation to observed difference in OA signal (ΔS) in Eq. (4), which can be substituted into Eq. (1) (See Eq. (5)). Simplification of the resulting equation, along with a major assumption on light fluence, yields an expression for HS680 kinetics in terms of OA signal rather than HS680 concentration (Eq. (6)). Light fluence between the aorta and tumor "region of interest" can be considered to be the same because both are approximately the same distance from the surface of the mouse, 2.51 mm and 2.35 mm (Fig. 3).

HS HS
680 680 HS680 was administered in the form of a 100 μL single bolus, IV injection via the tail vein catheter (100% bioavailability) of the female athymic mouse. Contrast signal in the plasma was measured in the aorta and quantified with the "region of interest" function in ViewMSOT software, and referred to as the arterial input function (AIF). Given the quantification of AIF, these known values were used to solve contrast concentration within the tumor microenvironment by substituting AIF for ΔS plasma in Eq. (6). A MATLAB function was written to solve the differential equation (Eq. (6)) for the contrast tumor time-concentration profile (C tumor ) given AIF OA signal intensity values and estimated rate constant parameters. The experimental MSOT data for HS680 within the tumor was used to estimate compartment model rate constant parameters, which was performed with the nonlinear least squares (lsqnonlin) algorithm in MATLAB. Lastly, the algorithm calculates the correlation value (R value) for the PK model against the experimental C tumor data.

Isolating target absorber in images
Biodistribution of HS680 within the breast cancer mouse model using MSOT imaging was visualized via multispectral unmixing. This process separated the unique HS680 spectral signature peaked at 680 nm from the background signals, identifying HS680 signal and differentiating signal intensity within the animal. Cross-sectional images of the mouse were taken during the imaging process at the location of the breast tumor. The first row showed MSOT images of grayscale background. The second row showed the intensity of oxygenated hemoglobin (HbO 2 ) biodistribution within the animal. The third row showed the intensity of HS680 biodistribution. The third row showed an image overlay of HS680 signal gradient on top of a singlewavelength background (Fig. 4). The background optoacoustic image was produced by a single-wavelength signal at 900 nm. HbO 2 reading provided the optimal method to define the anatomical structures of the animal. The optoacoustic image overlay process visualizes HS680 biodistribution and signal intensity with respect to live mouse anatomy.

Pharmacokinetic imaging for intratumoral HypoxiSense 680 dynamics
MSOT imaging was selected for continuous imaging of HS680 dynamic within the tumor microenvironment of a living mouse. For the purpose of the intratumoral PK study, cross-sectional MSOT images were captured in the region of the tumor over time after intravenous injection of HS680. Spectral unmixing of the images allowed for isolation of contrast signals. Overlaying these signals on top of the singlewavelength grayscale background provided visualization of HS680 kinetics within the tumor over time. Signal increased in the tumor and the tumor environment as time progressed. Accumulation of HS680 was confirmed with near infrared fluorescent imaging using Advanced Molecular Imager (AMI-1000X) 24 h post-injection ( Supplementary  Fig. 1).
"Region of interest" (ROI) tool in the ViewMSOT software was used to define regions. The three regions identified (Fig. 3) represented the tumor, aorta, and a region outside of the mouse. Defining ROIs allowed for calculation of the mean contrast signal intensity within these regions using spectral unmixing. The ROIs shown have areas 3.01 mm 2 for the tumor ROI and 0.82 mm 2 for the AIF ROI. "Region of interest" of the AIF and tumor were drawn with the intention to keep the area of the corresponding ROIs constant in other mice. The empirical data points calculated from the tumor ROI represented the experimental concentration-time profile for HS680 within the tumor of the mouse model. A 2-compartmental PK algorithm was written using MATLAB in which Fig. 2. The pharmacokinetic model schematics. The basic structure of the PK model consists of two compartments: plasma compartment and tumor compartment. The two parameters, k in and k out , represent the rates at which a contrast agent fluxes in and out of the tumor compartment, respectively. Based on the pictorial schematic, the two-compartmental PK model assumes homogeneous distribution and constant volume in each compartment. Parameters are derived from curve-fitting the PK model curve to empirical MSOT data.  the compartments were used to account for blood plasma and tumor microenvironment. The model was fitted to the empirical data points by optimizing the relevant rate constants (Fig. 5).

Discussion
The non-invasive nature and capability to perform live imaging with MSOT allow advanced novel tracking capabilities of intratumoral contrast dynamic and biological processes within a native environment. Decreasing variability between animals combined with the individual AIF model provides more potential patient-specific applications in the future. MSOT image acquisition is in the range of 0.1 s while MRI image acquisition is in the range of 10 s to 100 s [36,37]. Significantly faster image acquisition makes MSOT a better imaging option for a PK study by providing more data points in a given period of time (Fig. 5). In addition, the lower level of photon scattering from acoustic waves contributes to high spatial resolution at clinically relevant depth compared to optical imaging, computed tomography, and MRI [1].
MSOT images of the same slice of tumor at different time points after HS680 injection (0, 5, 10, 15, 20, 25 min.). It demonstrates the unique advantage of collecting information in regard to multiple chromophores acquired during a single scan [12]. The overlay provides spatial confirmation of the location of the tumor as the HS680 signal is located where HbO 2 signal is depleted. In addition, visualization of the increased and sustained HS680 signal at the tumor location confirms the quantified signal data obtained for the PK study. This study with HS680 provides previously unknown quantification of intratumoral kinetics with hypoxia-targeting agents, which could be compared with in vivo kinetics of hypoxic-targeting therapeutics and enable further understanding of dosing and time-concentration profiles of therapeutic agents.
Another unique aspect of this investigation is the utilization of HS680 to target carbonic anhydrase IX, which regulates hydrogen ion flux in cell membrane and is a known tumor hypoxia marker [15]. This investigation identifies the HS680 kinetics and accumulation within the tumor in comparison to the plasma (Fig. 3). Selection of a hypoxia targeting contrast agent differentiates this study from previous PK studies performed with nonspecific blood-pool contrast agents, which are more suitable for applications such as vascular imaging [38] and lymph node mapping [39].
Selection of HS680 as the contrast agent contributes to the time delay observed (Fig. 5). The increased vascular permeability based on the enhanced permeability and retention (EPR) effect in most solid tumors indicates that diffusion is the main method of contrast transport into the tumor microenvironment [40]. Thus, size is a main determinant of kinetics within the tumor. HS680 consists of a modified acetazolamide conjugated to IR-680 fluorochrome, having a molecular weight of 1500 g/mol. In comparison, ICG has a molecular weight of 774.96 g/mol while gadolinium-based MRI contrast agents have molecular weights ranging from 36 to 480 g/mol [41]. Thus, HS680 diffusion into the tumor microenvironment is delayed compared with these commonly used smaller contrast agents. This possible explanation is affirmed by previous PK studies using ICG and MSOT imaging, but revealed no time delay between tumor and AIF signals [18,26]. This 2-compartmental PK model requires only two rate parameters, contrast influx into (k in ) and efflux out (k out ) of tumor microenvironment (Fig. 2). The model in this investigation was constructed using a similar AIF approach as the Flow-Limited Model for DCE-MRI [33], which calculated HS680 concentration within the tumor microenvironment. Avoiding possible error derived from overfitting due to unnecessary parameters is the main reason for selecting a 2-compartmental PK model over a more detailed and complicated PBPK model. A PBPK model consists of multiple compartments that correspond to different tissues in the body, which are not necessary in this specific kinetic modeling approach using an AIF [42].
The AIF approach is utilized in quantitative DCE-MRI analyses. It uses an input measurement to describe arterial concentration of a contrast agent. An image-derived AIF approach must satisfy two criteria: visibility in target field of view and sufficiently high temporal resolution [43]. Alternative approaches have produced similar results using a population AIF in preclinical setting [44] and clinical setting [45,46] or using the reference region (RR) approach [47][48][49], which replaces input AIF with a differential equation that quantifies contrast agent accumulation in a well-characterized tissue (i.e. muscle) over time. The individual AIF approach is selected for this study with MSOT because it generates the most accurate data for each individual animal and it satisfies both crucial criteria for using individual AIF approach. First, the aorta was selected because the vessel was present in the same field of view as the tumor (Fig. 3). More importantly, MSOT has the temporal resolution required to capture the rapid dynamic of the contrast agent for PK curve characteristics. Both the population AIF or RR approach are alternative approaches in DCE-MRI because of limitations to MRI acquisition time. MSOT avoids this temporal resolution limitation and can utilize an individual AIF approach. The modeling results of the AIF approach are shown in Fig. 5. The deviation between PK model curve and experimental data is likely due to the high noise level of the AIF (Fig. 5). The AIF model quantified signal from the ROI representing the mouse aorta, which is prone to changes in vessel dimension, resulting in high level of noise. Breathing could have contributed to the cyclic noise patterns through flux of oxygenated and deoxygenated hemoglobin. Because the AIF is used directly to calculate the tumor concentration of HS680 (Eq. (1)), noise in the AIF is amplified in the PK model curve. Noise does not significantly affect the initial spike because of the signal spike from the bolus injection. A motion compensation algorithm can be applied to offset the rhythmic fluctuation of noise [50]. Another possibility is to attempt the reference region (RR) model rather than the individual AIF model. However, the RR method yielded significantly more error than the AIF method in a similar type of study using DCE-MRI [51].
The characterization of intratumoral kinetics using MSOT is not only helpful in tracking a contrast or drug in a mouse model but is also meaningful for clinical applications. Pharmacokinetic models have been developed to predict anticancer drug dosing and tumor growth [52], aide in development and optimization of anticancer therapy [53], and have been used to characterize antineoplastic drug interactions in human prostate cancer models [54]. We expect increasing clinical applications for the combination of target specific contrast agent and MSOT imaging such as for drug delivery. Numerous studies have investigated the effectiveness of various antineoplastic drug combinations in treating malignancy [55][56][57]. Drug regimens can benefit from the knowledge of peak intratumoral kinetics of the anticancer drugs, thus allowing for optimization of the synergistic effects of drug combinations. Another critical potential area of clinical application is for intraoperative detection of tumor for surgical resection [1]. Accurate prediction of peak intratumoral kinetics of a contrast agent during surgical resection would allow for improved surgical efficacy, minimize unnecessary surgical procedures, and increase the chances of successful tumor resection.

Conflict of interest
All contributors are associated with Wake Forest Baptist Comprehensive Cancer Center and Wake Forest University School of Biomedical Engineering.
Ted Xiao completed his Bachelor in Science at the California Institute of Technology and is a MD/PhD candidate at Wake Forest School of Medicine. He is currently a graduate student in the Virginia Tech-Wake Forest University School of Biomedical Engineering and Sciences. He is a member of the McNally Lab and focuses on nanoparticle construction and MSOT imaging in his ongoing PhD thesis. He is interested in effective cancer treatments and detection methods through application of nanotechnology and imaging techniques.   He is the author of over 200 peer-reviewed publications and has been awarded three patents on medical image analysis.
Umit Topaloglu, Ph.D., is the Associate Director for Informatics at the Wake Forest University Comprehensive Cancer Center and an Associate Professor of Cancer Biology and Biostatistics. He develops and implements a Biomedical Informatics Vision that supports and advances the clinical and translational research needs around the scalable, interoperable, findable, sharable data frameworks. He has worked on pathological and radiological image management and secure access including leading statewide systems for trauma and telemedicine network implementation. Additionally, he has focused automated or semi-automated identification of features utilizing variety of machine learning and deep learning methods in spatio-temporal domain, which has the promise for detecting agents to accomplish the proposed tasks. McNally's laboratory develops actively targeted contrast agents, theranostic nanoparticles, and novel treatments with the overall goal to improve patient outcomes. She uses multispectral optoacoustic tomography to track contrast agents and nanovehicles in both the pre-clinical and clinical setting.


extfieldsofstudy: Medicine Biology

text:
Complete plastome sequences of two Psidium species from the Galápagos Islands

We report the complete plastome sequences of an endemic and an unidentified species from the genus Psidium in the Galápagos Islands ( P. galapageium and Psidium sp. respectively).

Introduction
Over a quarter of all vascular plant species are endemic to islands, making them hotspots of plant diversity and conservation (Kreft et al., 2008). In the Galápagos Islands, there are roughly 560 native species of plants of which approximately 32% are endemic (Lawesson et al., 1987). However, many of these endemic species have remained relatively unstudied since they were originally given scientific descriptions, making the study of the evolutionary histories of these unique taxa difficult. In the present study, we constructed the complete plastome sequences of two species of Psidium (guava) from the Galápagos Islands, one endemic and one currently unidentified in hopes of facilitating future work on the evolutionary relationships of these species.

Methods
This research is authorized under the permit: MAE-DNB-CM-2016-004 in compliance with Ecuadorian regulations.
Leaf samples were collected during May of 2017 from the Galápagos endemic Psidium galapageium Hook (commonly known as guayabillo) on the island of San Cristobal (0.89094°S, 89.43769°W) and from an unidentified Psidium species on the island of Santa Cruz (0.62313°S, 90.38581°W). Based on morphological similarity, the Psidium sp. individual is suspected to be P. acidum (Landrum, 2016), but no reference or barcode sequence from P. acidum is available for confirmation.
Leaf tissue was desiccated immediately after harvesting using silica gel. DNA extractions were performed using a Qiagen DNeasy Plant mini kit (Qiagen, Inc.). Sequence data was generated in the form of paired-end, 150 bp reads using a KAPA library prep kit (Roche Sequencing) and sequenced on an Illumina HiSeq 4000 platform (Illumina, Inc.).
Reads were quality and adapter trimmed using Trim Galore! version 0.4.3 with a minimum phred score value of 20 and minimum read length of 50 bp. Filtered reads were then aligned to the Psidium guajava plastome reference available at NCBI (Accession: KX364403) using the mem function within BWA version 0.7.15 (Li & Durbin, 2009). Consensus plastome sequences were generated using the mpileup function within samtools version 1.8 followed by the call and consensus functions within bcftools with a minimum depth of coverage of 10x (Li et al., 2009). Using IRscope (Amiryousefi et al., 2018), the P. galapageium and Psidium sp. plastomes respectively were confirmed to contain a large single copy of 88,268 bp and 87,747 bp and a small single copy of 18,465 bp and 18,490 bp separated by two inverted repeats of 26,071 bp and 26,360 bp for total lengths of 158,875 bp and 158,957 bp (Figure 1).
Annotations were generated using the program Plann (Huang & Cronk, 2015). Of the 132 gene features annotated previously in the Psidium guajava (guava) chloroplast genome on NCBI (Accession: KX364403), all were recovered in the Psidium sp. and P. galapageium plastome sequences. The non-identity of the two taxa sampled is evidenced by the absolute pairwise sequence divergence of the concatenated sequences of three conserved genes (MatK, psbA, and rbcL) which have been successfully used as barcodes previously in Psidium (Kress et al., 2009). Sequences were aligned using MUSCLE within MEGA version 7.0.26 (Tamura et al., 2007), and the number of nucleotide differences were counted between these alignments to estimate divergence. A total of 35 differences were observed among 4011 sites (0.87% uncorrected divergence) between P. guajava (Accession: KX364403) and P. galapageium, 45 differences (1.1%) between P. guajava and Psidium sp., and 40 differences (0.99%) between P. galapageium and Psidium sp.

Data availability
Voucher specimens for P. galapageium and Psidium sp. are available at the Charles Darwin Research Station herbarium (Index Herbariorum code CDS) with accession numbers 3053515 and 3053562, respectively. The corresponding plastome sequences for P. galapageium and Psidium sp. are available at NCBI with accession numbers MH491846 and MH491847, respectively.

Competing interests
No competing interests were disclosed.

Grant information
This work was supported by a Louise Coker Fellowship from the University of North Carolina at Chapel Hill (UNC).
All funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


extfieldsofstudy: Medicine Chemistry

text:
Increased Hematopoietic Stem Cells/Hematopoietic Progenitor Cells Measured as Endogenous Spleen Colonies in Radiation-Induced Adaptive Response in Mice (Yonezawa Effect)

The existence of radiation-induced adaptive response (AR) was reported in varied biosystems. In mice, the first in vivo AR model was established using X-rays as both the priming and the challenge doses and rescue of bone marrow death as the end point. The underlying mechanism was due to the priming radiation-induced resistance in the blood-forming tissues. In a series of investigations, we further demonstrated the existence of AR using different types of ionizing radiation (IR) including low linear energy transfer (LET) X-rays and high LET heavy ion. In this article, we validated hematopoietic stem cells/hematopoietic progenitor cells (HSCs/HPCs) measured as endogenous colony-forming units-spleen (CFU-S) under AR inducible and uninducible conditions using combination of different types of IR. We confirmed the consistency of increased CFU-S number change with the AR inducible condition. These findings suggest that AR in mice induced by different types of IR would share at least in part a common underlying mechanism, the priming IR-induced resistance in the blood-forming tissues, which would lead to a protective effect on the HSCs/HPCs and play an important role in rescuing the animals from bone marrow death. These findings provide a new insight into the mechanistic study on AR in vivo.

Introduction
Ionizing radiation (IR) at high doses is detrimental to the exposed organism, and exposure to IR could increase the risk of developing cancers and other health-related problems even including acute radiation syndrome such as bone marrow death. However, biological effects of IR at low dose or low dose rate remain elusive. On the other hand, adaptation to the environmental genotoxic stresses or insults, such as IR, is one of the fundamental characteristics of life. In particular, prior mild stresses can provide some aid to prepare organisms for subsequent more severe stresses. 1 Radiation-induced adaptive response (AR) is a phenomenon that a priming low dose of IR induces resistance to a subsequent challenge exposure to IR at higher doses, and its studies could provide important scientific basis for IR risk estimates, protection, and practical applications. 2 Since the editio princeps of AR concept introduced into radiation biology, 3 it has been demonstrated in a variety of in vitro, in utero, and in vivo systems with end points such as DNA damage, chromosomal aberrations, cell transformation, cell death, and mutation in the in vitro experiments and prenatal death, malformation, hematopoietic death, and carcinogenesis. 1,[4][5][6] High atomic number and energy (HZE) particles such as carbon, oxygen, silicon, and iron are important component sof space radiation, from the solar particle events and the galactic cosmic rays. 7,8 Compared to photon and proton radiation, HZE particles bearing higher energy could cause both acute and long-term damage to bone marrow via increased production of reactive oxygen species, showing stronger detrimental effects (with higher relative biological effectiveness) on the hematopoietic system including decreased peripheral blood counts and reduced hematopoietic stem cells (HSCs) and progenitor cells (HPCs) in laboratory animal models [9][10][11][12][13][14][15][16][17][18][19][20] and, in addition, raising hematological cancer risk via bone marrow cell reprograming. 21 The AR mouse model for rescuing bone marrow death established by Yonezawa and colleagues 1,22-24 was repeatedly verified. 25,26 This model was named generally as "Yonezawa Effect" in Japan, which was originally established using low linear energy transfer (LET) X-rays as both the priming and the challenge IR, with the underlying mechanism that induction of radioresistance in blood forming tissues by the priming IR rescued the bone marrow death caused by the challenge high dose. In a series of investigations in our laboratory, we first verified and confirmed the existence of AR in mice using low-LET X-rays to deliver both the priming and the challenge IR and the 30-day survival test to estimate the efficacy for rescuing bone marrow death. Then, we further demonstrated the existence of AR using low-LET X-rays as the priming IR and high-LET heavy ion IR (carbon, neon, and silicon particles) as the challenge IR. Recently, we showed that, in this model, AR could be induced by high-LET heavy ion IR (carbon particles) as the priming IR and low-LET X-rays or high-LET heavy ion IR (carbon and neon particles) as the challenge IR. [20][21][22][23][24][25][26][27][28][29] In the present work, we validated the HSCs//HPCs measured as endogenous colony-forming units-spleen (CFU-S) under AR inducible and uninducible conditions using combination of low-LET X-rays and high-LET heavy ion IR as priming and challenge IR. By verifying the number of CFU-S under different conditions, 12 days after the animals received the challenge IR, the present investigation was aimed to study the recovery of HSCs/HPCs in this AR mouse model. We confirmed the consistency of CFU-S number with the AR induction conditions: significantly increased number of CFU-S under AR inducible conditions and no markedly increased number of CFU-S under AR uninducible conditions.

Animals
Five-week-old C57BL/6J Jms strain female mice were purchased from SLC, Inc (Japan) and maintained in a conventional animal facility under a 12 hour light-12 hour dark photoperiod. The animals were housed in autoclaved aluminum cages with sterilized wood chips and allowed to access standard laboratory chow (MB-1; Funabashi Farm Co, Japan) and acidified water ad libitum. The animals were acclimatized to the laboratory conditions for 1 week before use. To avoid possible effects from the developmental condition of the animals, 6-week-old mice with a significantly different body weight (more or less than the mean þ 2 SD) were omitted from this study. Based on our preliminary trials, in the present study, at least 6 mice were used in each experimental point. All experimental protocols involving mice were reviewed and approved by The Institutional Animal Care and Use Committee of the National Institute of Radiological Sciences (NIRS). The experiments were performed in strict accordance with the NIRS Guidelines for the Care and Use of Laboratory Animals.

Irradiation
For low-LET IR, X-rays were generated with an X-ray machine (Pantak-320S; Shimadzu, Japan) operated at 200 kVp and 20 mA, using a 0.50 mm Al þ 0.50 mm Cu filter. An exposure-rate meter (AE-1321M; Applied Engineering Inc, Japan) with an ionization chamber (C-110, 0.6 ml, JARP, Applied Engineering Inc, Japan) was used for the dosimetry. The dose rate for delivering the priming dose and the challenge dose was at about 0.30 Gy/min and 0.90 Gy/min, respectively. For high-LET heavy ion IR, the monoenergetic ion beam of carbon, neon, and iron particles was generated and accelerated by a synchrotron, the Heavy Ion Medical Accelerator in Chiba at NIRS, Japan. The beam energy was 290 MeV/nucleon, 400 MeV/nucleon, and 500 MeV/nucleon for carbon, neon, and iron particles, corresponding to an average LET value of about 15 keV/mm, 30 keV/mm, and 200 keV/mm, respectively. The dose rate was at about 0.10 Gy/min and 2.00 Gy/min for delivery of the priming dose and the challenge dose, respectively. The mice held in acryl containers were exposed to total body irradiation at room temperature.

Mouse Model for Radiation-Induced AR
The AR mouse model for rescue of bone marrow death and study on increase in the number of CFU-S established by Yonezawa and colleagues 24,30 was adopted, verified, and confirmed under the experimental conditions in our research facilities and finally applied to a series of our investigations using both low-LET X-rays and high-LET particles. The timing for delivery of the priming dose and challenge dose was on postnatal ages of 6 and 8 weeks of the mice, respectively. The present work was a part of the investigations focusing on the recovery of HSCs/HPCs under AR. The AR inducible and uninducible conditions obtained in our previous studies for rescue of bone marrow death [19][20][21][22][23][24][25][26][27] are summarized (Table 1). In Table 1, the "Yes" for AR induction meant a successful AR induction as judged by a significant increase in the animal survival after receiving the priming radiation prior to the challenge radiation in the 30-day survival test, and the "No" for no significant increase in the survival was induced in the presence of the priming radiation. To look for suitable experimental conditions, especially both the challenge dose altitude of each type of IR and the timing that makes the CFU-S more distinctly observable, 3 to 4 sublethal doses for the challenge IR of each type in each combination of the priming and the challenge exposure were validated in preliminary trials, and the doses listed in Table 1 were finally used. In brief, the combination exposures were (1) the priming dose was 0.50 Gy for low-LET X-rays and 0.45 Gy for high-LET carbon or iron particles; (2) the challenge dose for X-rays was 5.00 Gy in the combined exposure to priming X-rays and 4.75 Gy to priming carbon or iron particles; (3) the challenge dose for carbon particles was 5.00 Gy and 5.25 Gy following the priming IR from X-rays and carbon particles, respectively; and (4) the challenge dose was 5.00 Gy and 5.5 Gy for neon particles and iron particles, respectively.

Enumeration of Endogenous CFU-S
To evaluate the number of HSCs/HPCs in vivo in mice under experimental conditions capable or incapable of inducing AR, the method established by Till and McCulloch, 31 bearing advantages over the exogenous transplant system, including simplicity and rapidity, avoiding in vitro cell manipulation, 32,33 was adopted and applied to the present study. In brief, the mice were killed 11, 12, or 13 days after receiving the challenge dose. The spleens were removed, weighed, and then fixed in Bouin fixative (having a mixing ratio 15:5:1 of saturated picric acid, formaldehyde, and glacial acetic acid; purchased from Wako Pure Chemical Industries, Ltd, Japan) for 24 hours. The macroscopic nodules or colonies on the surface of the organ were counted as endogenous CFU-S using a stereoscopic microscope (Nikon SMZ-10; Nikon Instech Co, Ltd, Japan) at 10Â magnification. Timing of spleen sample collection was according to that the numeration of endogenous CFU-S formed 11 days or later after irradiation provides a measure of viable pluripotent HSCs/HPCs, 31,34 and longer intervals were avoided due to the colonies fused together, preventing scoring of individual colonies. Based on the ease for clear distinction of individual colonies, data obtained on 12 days after the challenge IR were used to judge the consistency of CFU-S number change and AR induction. Six to 12 animals were used per experimental point.

Statistical Analysis
Statistical evaluation of the data was done using the Student t test, and the statistical significance was assigned to P < .05.

Verification of the Radiation-Induced AR Mouse Model Using CFU-S as the End Point
Reproducibility of the radiation-induced AR mouse model (Yonezawa Effect) using CFU-S as the end point 30 was verified under the experimental setup in the present study. Under the AR inducible conditions, the animals were total body irradiated with a priming dose of 0.50 Gy X-rays at postnatal 6 weeks followed by a challenge dose of 5.00 Gy X-rays at postnatal 8¼weeks. Under the AR uninducible conditions, the animals were total body irradiated with only a challenge dose of 5.00 Gy X-rays at postnatal 8 weeks. The number of CFU-S was measured on 11, 12, and 13 days after the challenge IR. Results showed that the priming dose markedly increased the mean number of CFU-S from 3.2 to 11.8, 4.0 to 36.2, and 5.6 to 50.2 on the days 11, 12, and 13, respectively, after the challenge IR (Figure 1). Results clearly indicated that AR was induced with efficient reliability and reproducibility in our experimental setup using the number of CFU-S as the end point. Serving also as a positive control, the verification work was performed in parallel to the following investigations using combination of different types of IR.

Validation of CFU-S in AR Induced by Priming IR From Low-LET X-Rays and Challenge IR From High-LET Particles
The CFU-S assay was performed to validate whether significant increase in the number of CFU-S occurred in the animals under AR inducible condition (exposure of a priming dose of 0.50 Gy X-rays at postnatal 6 weeks followed by a challenge dose of 5.00 Gy carbon particles at postnatal 8 weeks) and the AR uninducible condition (the animals irradiated with only the challenge dose; the animals irradiated with a priming dose of 0.50 Gy X-rays at postnatal 6 weeks followed by a challenge dose of 5.50 Gy iron ions). The mean number of CFU-S was significantly increased from 7.8 to 15.0, 15.9 to 28.6, and 24.3 to 44.0, respectively, on 11, 12, and 13 days after the challenge IR ( Figure 2A). On the other hand, no increased number of CFU-S was observed in the animals receiving both the priming X-rays and the challenge iron IR when compared to the animals receiving only the challenge iron IR ( Figure 2B). On the 11th day after the challenge iron IR, the mean number of CFU-S was even markedly higher in the animals receiving only the challenge iron IR when compared to the animals receiving both the priming X-rays and the challenge iron IR ( Figure 2B), indicating an additive effect on reducing CFU-S from the combined exposure. These results clearly showed that under AR inducible condition when the priming IR was from low LET X-rays and the challenge IR was from high LET carbon particles, increased number of CFU-S was confirmed. On the other hand, under AR uninducible condition, when the priming IR was from low-LET X-rays and the challenge IR was from high-LET iron particles, no increased number of CFU-S was confirmed.  Figure 3A). The mean number of CFU-S was significantly increased from 6.7 to 22.3 and 4.2 to 32.5, respectively, on 12 and 13 days after the challenge IR ( Figure 3A). On the other hand, no increased mean number of CFU-S was observed in the animals receiving both the priming dose from iron particles followed by the challenge dose from X-rays and the animals irradiated with only the challenge dose from X-rays ( Figure 3B). These results clearly showed that under AR inducible condition when the priming IR was from high-LET carbon particles and the challenge IR was from low-LET X-rays, increased mean number of CFU-S was confirmed. Under AR uninducible condition, when the priming IR was from high LET iron particles and the challenge IR was from low LET X-rays, no increased mean number of CFU-S was observed.

Validation of CFU-S in AR Induced by High-LET Particles as Both Priming IR and Challenge IR
The CFU-S assay was performed to validate whether significant increase in the number of CFU-S occurred under AR inducible condition (ie, the animals irradiated with a priming dose of 0.45 Gy carbon particles at postnatal 6 weeks followed by a challenge dose of 5.25 Gy carbon particles or 5.00 Gy neon particles at postnatal 8 weeks vs the animals irradiated with only the challenge dose) and the AR uninducible condition (ie, the animals irradiated with a priming dose of 0.45 Gy carbon particles at postnatal 6 weeks followed by a challenge dose of 5.50 Gy iron particles at postnatal 8 weeks). On 11, 12, and 13 days after the challenge IR, the mean number of CFU-S was significantly increased from 2.6 to 9.2, 7.0 to 22.4, and 15.5 to 43.0 for the animals receiving the combined exposure to carbon particles when compared to those receiving only the challenge IR from carbon particles ( Figure 4A); the mean number of CFU-S was from 1.3 to 2.8, 4.0 to 8.6, and 4.7 to 10.1 for the animals exposed to the priming IR from carbon particles followed by the challenge IR from neon particles, and the increase was statistically significant on the 12th day ( Figure  4B). On the other hand, no increased mean number of CFU-S was observed in the animals receiving both the priming dose from carbon particles followed by the challenge dose from iron particles when compared to those receiving only the challenge dose from iron particles ( Figure 4C). These results clearly showed that under AR inducible condition when the priming The number of CFU-S was measured on the days 11, 12, and 13 after the challenge IR. Data of each experimental point were from 6 to 12 mice. Two asterisks (**) indicate statistically significant differences (P < .01) between the 2 groups that were compared. One asterisk (*) stands of statistically significant differences (P < .05) between the 2 groups that were compared. Two asterisks (**) indicate statistically significant differences (P < .01) between the 2 groups that were compared. C stands for carbon and Fe stands for iron. IR was from high-LET carbon particles and the challenge IR was from high LET carbon or neon particles, increased mean number of CFU-S was confirmed. Under AR uninducible condition, when the priming IR was from high LET carbon particles and the challenge IR was from high LET iron particles, no increased mean number of CFU-S was observed.

Discussion
A better understanding of AR and other nontargeted effects is needed to understand to which extent application of low-dose IR might be beneficial to humans. 35 To date, investigations using the AR mouse model (Yonezawa Effect) have obtained many substantial achievements in the study of radiationinduced AR at whole body level. In a series of comprehensive investigations, Yonezawa and colleagues verified the existence of AR under a variety of experimental conditions (ie, doses of priming and challenge IR, intervals between priming and challenge IR, and age and strain of the animals). 23 These efforts helped this AR mouse model to lay a cornerstone for in vivo AR research. Of note, based on the priming dose and the interval between priming and challenge exposures and the timing for delivery of the priming dose, 2 different phenotypes of AR were observed, involving different mechanisms; the first Data of each experimental point were from 6 to 12 mice. One asterisk (*) stands for statistically significant differences (P < .05) between the 2 groups that were compared. Two asterisks (**) indicate statistically significant differences (P < .01) between the 2 groups that were compared. C stands for carbon, Ne stands for neon, and Fe stands for iron.
phenotype was induced 2 weeks after a 0.3 to 0.5 Gy priming IR, which was due to Trp53-dependent radioresistance in blood-forming tissues, 26,36 and the second phenotype was observed 2 months after a 0.05 to 0.1 Gy priming exposure and the AR resulted from the interaction between blood-forming tissue and the central nervous system. 37,38 As rescue of bone marrow death is the basic criteria for judgment of a successful induction of AR in mice under Yonezawa Effect, studying the recovery of HSCs/HPCs is of critical significance from the point of view of mechanism research, the model for the first phenotype was applied to the present work to validate the HSCs/HPCs measured as endogenous CFU-S under AR inducible and uninducible conditions using combination of low-LET X-rays and high-LET heavy ion IR as priming and challenge IR. It is known that bone marrow, as the site in the body where self-renewal and differentiation of HSCs to mature blood cells mainly occurs, is extremely radiosensitive. Exposure to IR at high doses could devastate bone marrow leading to bone marrow death. In addition, hematopoietic capability is the most critical factor preventing radiation-induced bone marrow death. Even sublethal doses of IR could cause a decrease in hematopoietic cells and a deficit to bone marrow hematopoietic microenvironment. Ionizing radiation-induced decline in total bone marrow hematopoietic cells is accompanied with elevated adipocytes into the marrow cavity, leading consequently to the inhibition of bone marrow microenvironment recovery and hematopoiesis. 39 As the number of endogenous CFU-S could reflect both the number of HSCs/HPCs 32,40 and the environment for hematopoiesis, 41 endogenous CFU-S assay is capable of evaluating the hematopoietic capability.
In this AR mouse model, previous studies showed that successful induction of AR by priming low-LET X-rays or g-rays was regardless of the dose rate, 42 and mechanistic study suggested that priming IR-induced decreased p53, Bax, and apoptosis positive cell accumulation in the spleen might favor the recovery of hematopoietic function from challenge IR-induced acute injury, manifesting as stimulated recovery of spleen weight and endogenous CFU-S, contributing to a decrease in bone marrow death. 43 Studies on the protective effects induced by low doses of low-LET IR indicated that the underlying mechanisms included enhanced antioxidative capacities, increased cellular DNA repair capacity leading to such as reduced initial DNA damage in AR in mice in vivo, 44,45 and reduced cell death and mutations in vitro. 5,46 These induced responses were tightly conserved throughout evolution, being basic responses critical to life. 4 On the other hand, as high-LET IR from heavy particles induced biological effects qualitatively different from those induced by low-LET IR from photons, 47 for example, high-LET IR induced more clustered DNA damage and higher rates of residual chromatin breaks, 48 cellular radiosensitivity correlated with the frequency of residual chromatin breaks, 49 and the recovery ratio of the potentially lethal damage depended on the quality of IR. 50 In the present work, increased endogenous CFU-S was observed in the animals under AR inducible conditions, being consistently well with the 30-day survival results: Under AR inducible conditions for rescue of bone marrow death, AR manifested as significantly increased number of HSCs/HPCs could be induced regardless the priming low dose of IR was from low-LET X-rays or certain high-LET heavy ions such as carbon particles. These results indicate that induction of AR may protect the ability of animals to hinder the decline in the total HSCs/HPCs through possibly inducing radioresistance in the hematopoietic tissue and maintaining the hematopoietic microenvironment. These findings suggest that induction of AR by low-LET and certain high-LET heavy ions may share at least some mechanisms in common. In fact, mechanistic study in vitro in cultured human fibroblasts showed that gene expression profiles following low-LET g-rays and decays of high-LET like 125 I shared the majority of genes in common, indicating that both types of IR elicited similar signal transduction pathways, 51 and the extent of DNA clustered damage may not be the major factor modulating gene expression after exposure. 52 Low doses of low-LET X-rays were effective in reducing chromosomal aberrations and mutation frequency induced by high-LET IR. 2,53 These findings suggest that the biological defense mechanisms induced by prior low doses from either low-LET IR or high-LET IR may be considered as effective countermeasures, being sufficient enough against the damages caused by subsequent higher challenge doses from either low-LET or high-LET IR. On the other hand, it is noticed that heavy ions with higher atomic number and higher energy (ie, iron particles) failed to induce AR regardless of its use as priming or challenge IR. These findings also suggest that induction of AR may depend on the quality of IR at least to a certain extent. Is there a threshold for the atomic number and higher energy of the heavy ions to induce AR in this mouse model? More questions remain to be answered through further research on especially the underlying molecular mechanisms.
The priming dose used in the present work was higher than 100 mGy which was extensively used in the field of so-called low-dose research. It is known that in the experimental study on AR induction or hormesis, the low doses used in the in vivo systems are often relatively higher than that used in the in vitro systems, and the altitude of the dose seems to be dependent of the biosystems. On the other hand, when thinking about the clinical application of induction of AR or hormesis for the treatment of patients with cancer to protect the normal tissue from being damaged by radiotherapy at high doses, 0.50 Gy could be considered as a low dose. More importantly, application of AR or hormesis should be more practical and acceptable for most of the patients via research and development of medication based on the molecular mechanisms underlying induction of AR or hormesis rather than application via exposure of the patient to priming low dose of radiation.
Taken together, results (Table 2) in the present study showed that under AR inducible conditions, regardless of the quality of IR (ie, low LET and high LET, photons and particles) for the combination of the priming dose and the challenge dose, the priming IR could induce an increased number of HSCs/ HPCs as measured by the number of endogenous CFU-S, which may contribute to the rescue of bone marrow death. Results indicated the significance and possible application of AR to the reduction in acute radiation syndrome induced by high dose from either low-or high-LET IR. These findings bring new knowledge to the characterization of the Yonezawa Effect by providing a new insight into the mechanistic study on the hematopoietic system in the AR mouse model in vivo.


extfieldsofstudy:
Medicine
Computer Science

text:
GeneSeqToFamily: a Galaxy workflow to find gene families based on the Ensembl Compara GeneTrees pipeline

Abstract Background Gene duplication is a major factor contributing to evolutionary novelty, and the contraction or expansion of gene families has often been associated with morphological, physiological, and environmental adaptations. The study of homologous genes helps us to understand the evolution of gene families. It plays a vital role in finding ancestral gene duplication events as well as identifying genes that have diverged from a common ancestor under positive selection. There are various tools available, such as MSOAR, OrthoMCL, and HomoloGene, to identify gene families and visualize syntenic information between species, providing an overview of syntenic regions evolution at the family level. Unfortunately, none of them provide information about structural changes within genes, such as the conservation of ancestral exon boundaries among multiple genomes. The Ensembl GeneTrees computational pipeline generates gene trees based on coding sequences, provides details about exon conservation, and is used in the Ensembl Compara project to discover gene families. Findings A certain amount of expertise is required to configure and run the Ensembl Compara GeneTrees pipeline via command line. Therefore, we converted this pipeline into a Galaxy workflow, called GeneSeqToFamily, and provided additional functionality. This workflow uses existing tools from the Galaxy ToolShed, as well as providing additional wrappers and tools that are required to run the workflow. Conclusions GeneSeqToFamily represents the Ensembl GeneTrees pipeline as a set of interconnected Galaxy tools, so they can be run interactively within the Galaxy's user-friendly workflow environment while still providing the flexibility to tailor the analysis by changing configurations and tools if necessary. Additional tools allow users to subsequently visualize the gene families produced by the workflow, using the Aequatus.js interactive tool, which has been developed as part of the Aequatus software project.

Dr Nicola Soranzo
Abstract: Background Gene duplication is a major factor contributing to evolutionary novelty, and the contraction or expansion of gene families has often been associated with morphological, physiological and environmental adaptations. The study of homologous genes helps us to understand the evolution of gene families. It plays a vital role in finding ancestral gene duplication events as well as identifying genes that have diverged from a common ancestor under positive selection. There are various tools available, such as MSOAR, OrthoMCL and HomoloGene, to identify gene families and visualise syntenic information between species, providing an overview of syntenic regions evolution at the family level. Unfortunately, none of them provide information about structural changes within genes, such as the conservation of ancestral exon boundaries amongst multiple genomes. The Ensembl GeneTrees computational pipeline generates gene trees based on coding sequences and provides details about exon conservation, and is used in the Ensembl Compara project to discover gene families.

Findings
A certain amount of expertise is required to configure and run the Ensembl Compara GeneTrees pipeline via command line. Therefore, we have converted the command line Ensembl Compara GeneTrees pipeline into a Galaxy workflow, called GeneSeqToFamily, and provided additional functionality. This workflow uses existing tools from the Galaxy ToolShed, as well as providing additional wrappers and tools that are required to run the workflow.

Conclusions
GeneSeqToFamily represents the Ensembl Compara pipeline as a set of interconnected Galaxy tools, so they can be run interactively within the Galaxy's userfriendly workflow environment while still providing the flexibility to tailor the analysis by changing configurations and tools if necessary. Additional tools allow users to subsequently visualise the gene families produced by the workflow, using the Aequatus.js interactive tool, which has been developed as part of the Aequatus software project.

Abstract
Background Gene duplication is a major factor contributing to evolutionary novelty, and the contraction or expansion of gene families has often been associated with morphological, physiological and environmental adaptations. The study of homologous genes helps us to understand the evolution of gene families. It plays a vital role in finding ancestral gene duplication events as well as identifying genes that have diverged from a common ancestor under positive selection. There are various tools available, such as MSOAR, OrthoMCL and HomoloGene, to identify gene families and visualise syntenic information between species, providing an overview of syntenic regions evolution at the family level. Unfortunately, none of them provide information about structural changes within genes, such as the conservation of ancestral exon boundaries amongst multiple genomes. The Ensembl GeneTrees computational pipeline generates gene trees based on coding sequences and provides details about exon conservation, and is used in the Ensembl Compara project to discover gene families.

Findings
A certain amount of expertise is required to configure and run the Ensembl Compara GeneTrees pipeline via command line. Therefore, we have converted the command line Ensembl Compara GeneTrees pipeline into a Galaxy workflow, called GeneSeqToFamily, and provided additional functionality. This workflow uses existing tools from the Galaxy ToolShed, as well as providing additional wrappers and tools that are required to run the workflow.

Introduction
The phylogenetic information inferred from the study of homologous genes helps us to understand the evolution of gene families, which plays a vital role in finding ancestral gene duplication events as well as identifying regions under positive selection within species [1]. In order to investigate these low-level comparisons between gene families, the Ensembl Compara GeneTrees gene orthology and paralogy prediction software suite [2] was developed as a pipeline that uses TreeBest [3] [4] (part of TreeFam [5]) to find internal structural-level synteny for homologous genes. TreeBeST implements multiple independent phylogenetic methods and can merge their results in a consensus tree whilst trying to minimise duplications and deletions relative to a known species tree. This allows TreeBeST to take advantage of the fact that DNAbased methods are often more accurate for closely related parts of trees, while protein-based trees are better at longer distances.
The Ensembl GeneTrees pipeline comprises seven basic steps, starting from a set of protein sequences and performing similarity searching and multiple large-scale alignments to infer homology among them, using various tools: BLAST [6], hcluster_sg [7], T-Coffee [8], and phylogenetic tree construction tools, including TreeBeST. Whilst all these tools are freely available, most are specific to certain computing environments, are only usable via the command line, and require many dependencies to be fulfilled. Therefore, users are not always sufficiently expert in system administration in order to install, run, and debug the various tools at each stage in a chain of processes. To help ease the complexity of running the GeneTrees pipeline, we have employed the Galaxy bioinformatics analysis platform to relieve the burden of managing these system-level challenges.
Galaxy is an open-source framework for running a broad collection of bioinformatics tools via a user-friendly web interface [9]. No client software is required other than a recent web browser, and users are able to run tools singly or aggregated into interconnected pipelines, called workflows. Galaxy enables users to not only create, but also share workflows with the community. In this way, it helps users who have little or no bioinformatics expertise to run potentially complex pipelines in order to analyse their own data and interrogate results within a single online platform. Furthermore, pipelines can be published in a scientific paper or in a repository such as myExperiment [10] to encourage transparency and reproducibility.
In addition to analytical tools, Galaxy also contains plugins [11] for data visualisation. Galaxy visualisation plugins may be interactive and can be configured to visualise various data types, for example, bar plots, scatter plots, and phylogenetic trees. It is also possible to develop custom visualisation plugins and easily integrate them into Galaxy. As the output of the GeneSeqToFamily workflow is not conducive to human readability, we also provide a data-tovisualisation plugin based on the Aequatus software [12]. Aequatus.js [13] is a new JavaScript library for the visualisation of homologous genes, which is extracted from the standalone Aequatus tool. It provides a detailed view of gene structure across gene families, including shared exon information within gene families alongside gene tree representations. It also shows details about the type of interrelation event that gave rise to the family, such as speciation, duplication, and gene splits.

Methods
The GeneSeqToFamily workflow has been developed to run the Ensembl Compara software suite within the Galaxy environment, combining various tools alongside preconfigured parameters obtained from the Ensembl Compara pipeline to produce gene trees. Among the tools used in GeneSeqToFamily (listed in Table 1), some were existing tools in the Galaxy ToolShed [14], such as NCBI BLAST, TranSeq, Tranalign and various format converters. Additional tools that are part of the pipeline were developed at the Earlham Institute (EI) and submitted to the ToolShed, i.e. BLAST parser, hcluster_sg, hcluster_sg parser, T-Coffee, TreeBeST best and Gene Alignment and Family Aggregator. Finally, we developed helper tools that are not part of the workflow itself, but aid the generation of input data for the workflow and these are also in the ToolShed, i.e. Get features by Ensembl ID, Get sequences by Ensembl ID, Select longest CDS per gene, ETE species tree generator and GeneSeqToFamily preparation.
The workflow comprises 7 main steps, starting with translation from input coding sequences (CDS) to protein sequences, finding subsequent pairwise alignments of those protein sequences using BLASTP, and then the generation of clusters from the alignments using hcluster_sg. The workflow then splits into two simultaneous paths, whereby in one path it performs the multiple sequence alignment (MSA) for each cluster using T-Coffee, and in the other, generates a gene tree with TreeBeST taking the cluster alignment and a species tree as input. Finally, these paths merge to aggregate the MSA, the gene tree and the gene feature information (transcripts, exons, and so on) into an SQLite [15] database for visualisation and downstream reuse. Each step of the workflow along with the data preparation steps is explained in detail below.   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 Get features by Ensembl ID

Data generation and preparation
We have developed a number of tools that assist in preparing the datasets needed by the workflows.
Ensembl REST API tools Galaxy tools were developed which utilise the Ensembl REST API [31] to retrieve sequence information (Get sequences by Ensembl ID) and feature information (Get features by Ensembl ID) by Ensembl ID from the Ensembl service. REST (REpresentational State Transfer) is an architecture style for designing networked applications [32] which encourages the use of standardised HTTP technology to send and receive data between computers. As such, these tools are designed to help users to retrieve existing data from Ensembl rather than requiring 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 them to manually download datasets to their own computers and then subsequently uploading them into the workflow.
We have also developed: • Select longest CDS per gene, which filters a CDS FASTA file from Ensembl retaining only the longest CDS sequence for each gene • ETE species tree generator, which uses the ETE toolkit [33] to generate a species tree from a list of species names or taxon IDs through NCBI Taxonomy.

GeneSeqToFamily preparation
Before GeneSeqToFamily can be run, a data preparation step must be carried out. We have developed a tool called GeneSeqToFamily preparation to preprocess the input datasets (gene feature information and CDS) for the GeneSeqToFamily workflow. It converts a set of gene feature information files in GFF3 [34] and/or JavaScript Object Notation (JSON) [35] format to an SQLite database. It also modifies all CDS FASTA header lines by appending the species name to the transcript identifier, as required by TreeBeST best. We decided to use an SQLite database to store the gene feature information because: • the GFF3 format has a relatively inconvenient and unstructured additional information field (9 th column) • searching is much faster and more memory efficient in a database than in a text file like JSON or GFF3, especially when dealing with feature information for multiple large genomes

CDS translation
Transeq Transeq, part of the European Molecular Biology Open Software Suite (EMBOSS) [36], is a tool to generate six-frame translation of nucleic acid sequences to their corresponding peptide sequences. Here we use Transeq to convert a CDS to protein sequences in order to run BLASTP and find protein clusters. However, since downstream tools in the pipeline such as TreeBeST require nucleotide sequences to generate a gene tree, the protein sequences cannot be directly used as workflow input and are instead generated with Transeq.

Pre-clustering alignment
BLAST This workflow uses the BLAST wrappers [37] developed to run BLAST+ tools within Galaxy. BLASTP is run over the set of sequences against the database of the same input, as is the case with BLAST-all, in order to form clusters of related sequences.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 BLAST parser BLAST parser is a small Galaxy tool to convert the BLAST output into the input format required by hcluster_sg. It takes the BLAST 12-column output [38] as input and generates a 3-column tabular file, comprising the BLAST query, the hit result, and the edge weight. The weight value is simply calculated as minus log10 of the BLAST e-value, replacing this with 100 if this value is greater than 100. It also removes the self-matching BLAST results.

Cluster generation
hcluster_sg hcluster_sg performs hierarchical clustering under mean distance for sparse graphs. It reads an input file that describes the similarity between two sequences, and iterates through the process of grouping two nearest nodes at each iteration. hcluster_sg outputs a single list of gene clusters, each comprising a set of sequence IDs present in that cluster. This list needs to be reformatted using the hcluster_sg parser tool in order to be suitable for input into T-Coffee and TreeBeST (see below).
hcluster_sg parser hcluster_sg parser converts the hcluster_sg output into a collection of lists of IDs, one list for each cluster. Each of these clusters will then be used to generate a gene tree via TreeBeST. The tool can also filter out clusters with a number of elements outside a specified range. The IDs contained in all discarded clusters are collected in separate output dataset. Since TreeBeST requires at least 3 genes to generate a gene tree, we configured the tool to filter out clusters with less than 3 genes.
Filter by FASTA IDs, which is available from the Galaxy ToolShed, is used to create separate FASTA files using the sequence IDs listed in each gene cluster.

Cluster alignment
T-Coffee T-Coffee is a MSA package, but can also be used to combine the output of other alignment methods (Clustal, MAFFT, Probcons, MUSCLE) into a single alignment. T-Coffee can align both nucleotide and protein sequences [8], and we use it to align the protein sequences in each cluster generated by hcluster_sg.

Gene tree construction
Tranalign Tranalign [36] is a tool that reads a set of nucleotide sequences and a corresponding aligned set of protein sequences and returns a set of aligned nucleotide sequences. Here we use it to generate CDS alignments of gene sequences using the protein alignments produced by T-Coffee.
TreeBeST 'best' TreeBeST (Tree Building guided by Species Tree) is a tool to generate, manipulate, and display phylogenetic trees and can be used to build gene trees based on a known species tree.
The 'best' command of TreeBeST builds 5 different gene trees from a FASTA alignment file using different phylogenetic algorithms, then merges them into a single consensus tree using a species tree as a reference. In GeneSeqToFamily, TreeBeST best uses the nucleotide MSAs generated by Tranalign (at least 3 sequences are required) and a user-supplied species tree in Newick format [39] (either produced by a third-party software or through the ETE species tree generator data preparation tool, described above) to produce a GeneTree for each family represented also in Newick format. The resulting GeneTree also includes useful annotations specifying phylogenetic information of events responsible for the presence/absence of genes, for example, 'S' means speciation event, 'D' means duplication, and 'DCS' denotes the duplication score.

Aequatus visualisation plugin
The SQLite database generated by the GAFA tool can be rendered using a new visualisation plugin, Aequatus.js. The Aequatus.js library, developed at EI as part of the Aequatus project, has been configured to be used within Galaxy to visualise homologous gene structure and gene family relationships. This allows users to interrogate not only the evolutionary history of the gene family but also the structural variation (exon gain/loss) within genes across the phylogeny. Aequatus.js is available to download from GitHub [41], as visualisation plugins cannot yet be submitted to the Galaxy ToolShed.

Finding homology information for orphan genes
Although the GeneSeqToFamily workflow will assign most of the genes to orthogroups, many genes within a species might appear to be unique without homologous relationship to any other genes from other species. This observation could be the consequence of the parameters selected, choice of species, incomplete annotations. This could also reflect real absence of homology such as for rapidly evolving gene families. In addition to the GeneSeqToFamily workflow, we also developed two associated workflows to further annotate these genes by: 1) Retrieving a list of orphan genes from the GeneSeqToFamily workflow (see Figure 5) as follows: a) Find the IDs of the sequences present in the input CDS of the GeneSeqToFamily workflow, but not in the result of BLAST parser from the same workflow b) Add to this list the IDs of the sequences discarded by hcluster_sg parser c) From the input CDS dataset, retrieve the respective sequence for each CDS ID (from the step above) using Filter by FASTA IDs These unique CDS can be fed into the SwissProt workflow below to find homologous genes in other species. 2) Finding homologous genes for some genes of interest using SwissProt (see Figure 6) as follows: a) Translate CDS into protein sequences using Transeq b) Run BLASTP for the protein sequences against the SwissProt database (from NCBI) c) Extract UniProt IDs from these BLASTP results, using the preinstalled Galaxy tool Cut columns from a table (tool id Cut1) d) Retrieve Ensembl IDs (representing genes and/or transcripts) for each UniProt ID using UniProt ID mapping and retrieval e) Get genomic information for each gene ID and CDS for each transcript ID from the core Ensembl database using Get features by Ensembl ID and Get sequences by Ensembl ID respectively. The results from this second workflow can be subsequently used as input to GeneSeqToFamily for familial analysis.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 Example use case Since BLASTP plays a crucial role in determining gene families, we tested this workflow on a large dataset of CDS from three vertebrate species in order to set a benchmark and find optimum parameters to run the workflow. We downloaded the CDS sequences for Sarcophilus harrisii (Tasmanian devil), Mus musculus (Mouse), and Ornithorhynchus anatinus (Platypus) from Ensembl (release 87) and filtered them to retain only the longest transcript per gene (as in the Ensembl Compara pipeline), obtaining a total of 62,597 CDS. We then ran the GeneSeqToFamily workflow on them using various BLASTP parameters (as shown in Table 3), in order to identify the optimal values for the workflow to generate a gene tree in which members are possibly evolved from a single ancestor gene, and usually with identical biochemical functions such as proteins.
Our results show that the number of gene families can vary quite distinctly with different BLASTP parameters. Stringent parameters (Analysis 6) result in a large number of smaller families, while relaxed parameters (Analysis 1) generate a smaller number of large families, which may include distantly related genes. By testing different parameters and comparing the analyses with third party tools such as PantherDB to validate the results against known families, we chose those parameters listed as Analysis 5. These values seem to consistently generate legitimate sets of gene families with closely related family members based on the datasets we tested.
There are caveats, however. BLASTP parameters used in Analysis 5 restrict the maximum number of target sequences per query (max_target_seqs) to 3 (the first hit when using allversus-all BLAST will always be the query sequence itself). The minimum query coverage per High-scoring Segment Pair (HSP) (qcovhsp) is set to 90% and e-value cut-off to 1e-10, in order to find the HSP closest to the query thus allowing partial matches which could be seen in the event of gene split. If the input CDS contain multiple alternative transcripts per gene, we recommend setting the max_target_seqs parameter to 4 rather than 3 to get a wider range of results from BLAST, thereby helping to generate gene families with matching genes together with alternative transcripts. In contrast setting a value of 3 for the max_target_seqs parameter will restrict the search to only 3 matches per query, and the presence of alternative transcripts will decrease the likelihood of finding matches from other genes, thus increasing the likelihood of splitting of a gene tree into multiple trees and adversely inflating the number of families. Table 2: Results of the GeneSeqToFamily workflow run on 62,597 CDS from 3 species using 6 different BLAST parameter configurations, the complete list of which are shown in Table 3.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63   Before running the workflow, feature information and CDS for the selected genes were retrieved from the core Ensembl database using the helper tools described above (Get features by Ensembl ID and Get sequences by Ensembl ID respectively). A species tree was generated using ETE species tree generator and CDS were prepared with GeneSeqToFamily preparation.
We ran the GeneSeqToFamily workflow on these data using the parameters shown in Analysis 5 of Table 3, but we set max_target_seqs as 4 (as described in the previous use case) to get a wider range of results from BLAST because our dataset includes alternative transcripts. This workflow generated 4 different gene trees, one for each gene family. Figure 7,

Conclusion
The ultimate goal of the GeneSeqToFamily is to provide a user-friendly workflow to analyse and discover homologous genes using the Ensembl Compara GeneTrees pipeline within the Galaxy framework, where users can interrogate genes of interest without using the command-line whilst still providing the flexibility to tailor analysis by changing configurations and tools if necessary. We have shown it to be an accurate, robust, and reusable method to elucidate and analyse potentially large numbers of gene families in a range of model and non-model organisms. The workflow stores the resulting gene families into a SQLite database, which can be visualised using the Aequatus.js interactive tool, as well as shared as a complete reproducible container for potentially large gene family datasets.
Gradually, we hope that the Galaxy community will undertake their own analyses and feedback improvements to various tools, and publish successful combinations of parameters used in the    GeneSeqToFamily workflow. We encourage this process by allowing users to share their own version of GeneSeqToFamily workflow for appraisal by the community.

Future directions
In terms of core workflow functionality, we would like to incorporate pairwise alignment between pairs of genes for closely related species in addition of the MSA for the gene family, which will help users to compare orthologs and paralogs in greater detail.
We also plan to include explicit integration of the PantherDB resources [47]. Protein ANalysis Through Evolutionary Relationships (PANTHER) is a classification system to characterise known proteins and genes in order to certify genomic annotation. Association of PantherDB with GeneSeqToFamily will enable the automation of gene family validation and add supplementary information about those gene families, which could then be used in turn to further validate novel genomics annotation.
We also plan to add the ability to query the GAFA SQLite database using keywords, to make it easy for users to find gene trees which include their genes of interest without needing to delve into the database itself.