DOI: 10.5061/dryad.brv15dv5x
title: DNA methylation profiling and genomic analysis in 20 children with short stature who were born small-for-gestational age

description: Purpose: In a significant proportion of children born
 small-for-gestational age (SGA) with failure of catch-up growth, the
 etiology of short stature remains unclear after routine diagnostic
 work-up. We wanted to investigate if extensive analysis of the (epi)genome
 can unravel the cause of growth failure in a significant portion of these
 children. Patients and Methods: Twenty SGA children treated with growth
 hormone (GH) because of short stature were selected from the BELGROW
 database of the Belgian Society for Pediatric Endocrinology and
 Diabetology for exome sequencing, SNP array and genome-wide methylation
 analysis to identify the (epi)genetic cause. First year response to GH was
 compared to the response of SGA patients in the KIGS database. Results: We
 identified (likely) pathogenic variants in 4 children (from 3 families)
 using exome sequencing and found pathogenic CNV in 2 probands using SNP
 array. In a child harboring a NSD1-containing microduplication, we
 identified a DNA methylation signature that is opposite to the genome-wide
 DNA methylation signature of Sotos syndrome. Moreover, we observed
 multi-locus imprinting disturbances in two children in whom no other
 genomic alteration could be identified. Five out of 6 children with a
 genetic diagnosis had an "above average" response to GH.
 Conclusions: The study indicates that a more advanced approach with deep
 genotyping can unravel unexpected (epi)genomic alterations in SGA children
 with persistent growth failure. Most SGA children with a genetic diagnosis
 had a good response to GH treatment.
description: <i>Patient recruitment</i> BELGROW, the Belgian database of children treated with growth hormone (GH), was searched for children born SGA and with short stature at the start of GH treatment. Other inclusion criteria were: (1) both biological parents are alive (2) bone age at the start of GH treatment is within 1.5 year of the calendar age; (3) Height SDS at the start of GH treatment is at least 1.5 SD below the height SDS of each parent; (4) Still prepubertal during the first year of GH treatment to allow a reliable assessment of the first year growth response to GH administration. Exclusion criteria were: (1) a known bone dysplasia or sitting height/total height &gt; 2SDS; (2) evidence for fetotoxic factors that explained the intra-uterine growth retardation (e.g. tobacco or alcohol abuse) or another known factor that explained the growth deficit or that modulated the response to GH therapy such as a chronic disease or chronic medication. Birth weight and length were calculated using the Niklasson references<sup>13</sup>. Height and weight SDS were calculated using the 2004 Flemish growth references<sup>14</sup>. Out of 65 eligible patients, 20 participated in the study. Reasons for non-participation included: one parent not available (died, in prison, moved to another country, …), patient had moved and address was unknown, incorrect registry data that when corrected turned the patient no longer eligible, genetic diagnosis reached but not available in the registry, refusal to participate, patient did not start treatment, patient was not invited for the study. <i>Ethical and regulatory aspects</i> The study was approved by the Academic Ethical Committee of the Brussels Alliance for Research and Higher Education (B200-2014-043). Both parents gave their informed consent and each child received an age-appropriate study information document and signed an assent form. It was stipulated in the study information document that secondary findings would not be communicated to the parents unless the genetic finding was featured on the list of the American College of Genetics and Genomics (ACMG) of genes, conditions and variants that are recommended to be reported back because of their important consequences for childhood health<sup>15</sup>. Alternatively, parents could opt out to be informed at all. <i>Evaluation of response to growth hormone</i> Height data at the start of treatment and 9-15 months later were extracted from the registry and scaled by intra- or extrapolation to 12 months. The height velocity during the first year of treatment was compared with the published height velocity response curves from a large registry (KIGS)<sup>16</sup>. The growth response of the patients was categorized as follows: “non-responder”: height velocity &lt; -1 SDS, “below average”: between -1 and 0 SDS, “above average”: between 0 and +1 SDS and “super-responder”: above +1 SDS. <i>Exome sequencing</i> DNA from 20 trios of children and their biological parents was extracted from peripheral white blood cells using standard methods. Sample preparation was done using the SeqCap EZ Library SR protocol (User’s Guide version 5.0, Roche). Paired-end sequencing (2x100bp or 2x150bp) was performed on the Hiseq 1500<sup>17</sup> (Illumina, San Diego, CA, USA) or the Nextseq 500<sup>18</sup> (Illumina, San Diego, CA, USA) sequencing platform based on availability. Data-analysis was performed using an in-house developed automated Galaxy based pipeline<sup>19</sup>. Variant calling was done using the Genome Analysis Toolkit Unified Genotyper<sup>20</sup><sup>,</sup><sup>21</sup>. Variants were annotated using ANNOVAR<sup>22</sup><sup>,</sup><sup>23</sup> in VariantDB<sup>24</sup><sup>,</sup><sup>25</sup>. GRCh37 was used as a reference build. <i>Automated variant prioritization </i> Automated variant prioritization was performed with the MOON interpretation software<sup>26</sup> (version 2.1.3, Diploid Orbicule BVBA, Heverlee, Belgium), as described previously<sup>27</sup>. In short, MOON uses artificial intelligence to prioritize genetic variants based on patient input data and a disorder model that consists of associations between diseases, disease genes and inheritance patterns. The disorder model is created and updated on a regular basis by performing natural language processing of the medical literature<sup>27</sup>. The following Human Phenotype Ontology (HPO) terms were used as input: short stature and/or proportionate short stature and/or small for gestational age and/or intrauterine growth retardation. Depending on the presence of additional clinical features, more patient-specific HPO terms were included. MOON listed variants showing autosomal dominant, autosomal recessive, X-linked dominant, or X-linked recessive inheritance that co-segregated with the phenotype in the family. Also reduced penetrance was taken into account. Variants were manually inspected with the Integrative Genomics Viewer<sup>28</sup><sup>,</sup><sup>29</sup> (IGV) to exclude false positives. Variant classification was performed using the criteria of the American College of Medical Genetics and Genomics (ACMG) and the Association for Molecular Pathology<sup>30</sup>. Variants were considered causal if the ACMG criteria predicted the variant to be pathogenic or likely pathogenic and if the identified variant was co-segregating with the phenotype in the family. Where necessary, this analysis was followed by a phenotypic re-evaluation of the family. <i>Variant validation by PCR and Sanger sequencing</i> Primers were designed using Primer3<sup>31</sup><sup>,</sup><sup>32</sup> (version 4.1.0). PCRs were performed using the GoTaq G2 kit (Promega, Leiden, the Netherlands) and checked on a 1 % agarose gel. PCR cleanup was performed using alkaline phosphatase (Roche, Basel, Switzerland) and exonuclease I (Bioké, Leiden, the Netherlands). Sanger sequencing was performed using the ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems Inc., Foster City, USA) and the ABI3130XL sequencer<sup>33</sup> (Applied Biosystems Inc., Foster City, USA). CLC DNA Workbench 5.0.2<sup>34</sup> (CLC bio, Aarhus, Denmark) was used for data analysis. <i>SNP array</i> Single nucleotide polymorphism (SNP) array analysis of the 20 children and their parents was performed using a HumanCytoSNP-12 v2.1 beadchip on an iScan system<sup>35</sup> (Illumina, San Diego, CA, USA), according to manufacturer instructions. Analysis of copy-number variations (CNV) was performed in the CNV webstore<sup>36</sup><sup>,</sup><sup>37</sup>. GRCh37 was used as a reference build. Homozygosity mapping was performed using Plink<sup>38</sup> in the CNV webstore. Classification of CNV was conducted according to the technical standards for the interpretation and reporting of CNV of the American College of Medical Genetics and Genomics (ACMG) and the Clinical Genome Resource (ClinGen)<sup>39</sup>. Where necessary, this analysis was followed by a phenotypic re-evaluation of the family. <i>Genome-wide methylation array </i> The Infinium MethylationEPIC BeadChip Kit (Illumina, San Diego, CA, USA) was used to target more than 850k CpGs in promoters, gene bodies, enhancers and intergenic regions in the genome (reference build GRCh37). DNA was bisulfite converted using the EZ DNA methylation kit (Zymo research, Irvine, CA, USA) according to manufacturer’s instructions. Bisulfite converted DNA was amplified, fragmented and hybridized onto the methylation BeadChips. After hybridization, unbound fragments were washed and hybridized fragments were extended using fluorescent nucleotide bases. For each probe present on the BeadChip, raw methylation data were obtained using the Illumina iScan system<sup>35</sup> (Illumina, San Diego, CA, USA) according to manufacturer instructions. For each CpG, the mean beta-value was calculated, where a beta-value of 0 represents no methylation and a beta-value of 1 represents full methylation. Details about the DNA methylation data analysis can be found in the supplementary methods<sup>40</sup>.