Michal Kraus

Lucia Šimová

Pavel Neogrády

Miroslav Urban
email: urban@fns.uniba.sk

Core -valence correlation effects in the ground and lowlying excited states of GaN

Keywords: GaN molecular properties, excited states, core correlation effects, CCSD(T), CASPT2, aug GaN molecular properties, excited states, core correlation effects, CCSD(T), CASPT2, aug-cc-pVXZ and aug-cc-PCVXZ basis sets

We present careful analysis of the accuracy of molecular properties of the GaN molecule calculated by the Coupled Cluster CCSD(T) method including the core-valence correlation effects. We demonstrate that calculations of the core-valence correlation effects with basis sets which were optimized considering just valence electrons of the gallium atom lead to unreliable results. Our objective is finding the relation between the accuracy and the oversimplification due to the use of aug-cc-pVXZ basis sets when 3d and eventually also 3s3p electrons are correlated with these basis sets, having in mind calculations of larger systems. We observe the oscillatory behaviour of molecular properties with increasing cardinal basis set number X in the aug-cc-pVXZ series. This is main deficiency of this more economical basis set. Extrapolation to the complete basis set limit with this series is only possible starting with the aug-cc-pVQZ basis. Regular behaviour with increasing X show results with the aug-cc-pCVXZ series. Reasonable molecular properties are provided by the ANO-RCC basis sets with 3d 10 4s 2 4p 1 correlated electrons of Ga.

Introduction

Saying that for reasonably accurate calculations of molecules containing a IIIA element the (n -1)d electrons must be correlated, along with three ns 2 np 1 valence electrons, means nothing new [1][2][3][4]. Less frequently is stressed that this should be accompanied by choosing a basis set, which is constructed having in mind correlation of deeper shells of a IIIA element [1,2]. In our recent paper [5] we have correlated at least (n-1)d 10 ns 2 np 1 electrons of gallium (or indium) and five 2s 2 2p 3 valence electrons of the nitrogen atom in calculations of molecular properties of the GaN and InN dimers. With the common aug-cc-pVXZ basis set series (X=T, Q, or 5) [6] we have observed unexpectedly large basis set superposition error (BSSE) [7] for some reaction energies corresponding to the breakdown of the (GaN) 2 or (InN) 2 dimer. This is a clear signal of the inadequacy of the basis set used for such calculations. Similar indication of errors due to the use of inappropriate basis set was carefully analyzed in calculations of GaN and other strongly bound systems, like GaNGa, by Tzeli and Tsekouras [2]. For molecules containing Ga to Kr elements DeYonker, Peterson and Wilson [1] have constructed and optimized the family of core -valence aug-cc-pCVXZ sets to correlate the 3s3p3d shell along with the valence 4s4p electrons within the spin-free Douglas-Kroll-Hess (DKH) relativistic approximation [8]. Such basis sets should be used for the treatment of the IIIA element containing molecules when correlating more than just valence ns 2 np 1 electrons of a metal element.

Clearly, aug-cc-pCVXZ basis sets are much larger than their precursors, the correlation consistent polarized valence aug-cc-pVXZ basis set family created by Wilson et.al. [6]. Combining large basis sets and large number of correlated electrons becomes prohibitively demanding when Coupled Cluster CCSD(T) or other highly sophisticated methods are used for calculations of compounds containing several Ga or other IIIA group elements. Therefore, for pragmatic reasons, it is tempting to correlate (n -1)d 10 ns 2 np 1 or more core-valence electrons, yet using a simpler aug-cc-pVXZ basis. However, aug-cc-pVXZ basis sets for Ga were constructed considering the correlation of the 4s 2 4p 1 valence shell only. Strictly speaking, correlating more than 3 electrons of a IIIA group element with these basis sets is inappropriate. Intuitively, one may feel that pushing such a basis sets "to do a job" for which they were not constructed, that is, to correlate (n -1)d 10 ns 2 np 1 or more electrons, is better than nothing and the error is hopefully small. DeYonker et. al [1] have concluded their investigation of this problem saying, "the detrimental effects of correlating the 3s3p3d electrons with just the valence aug-cc-pVXZ basis sets are clearly observed". Their observation is based on calculations of spectroscopic properties of the Ga 2 up to Kr 2 molecules in their ground state. Combining this warning with the fact that at least 3d electrons must be correlated for accurate (or even semi-accurate results) leaves a small space for a meaningful use of aug-cc-pVXZ basis sets for molecules containing Ga and similar elements. Unfortunatelly, the necessity of using aug-cc-pCVXZ basis sets may be quite demanding for large scale calculations.

Therefore, we have decided to perform CCSD(T) calculations of properties of the GaN molecule addressing two issues in this paper. The first one is a systematic investigation of the relation between the accuracy and the oversimplification due to using aug-cc-pVXZ basis sets when 3d and eventually also 3s3p electrons are correlated with these basis sets, having in mind calculations of larger systems. Along with careful consideration of the relation between the way in which the basis set was constructed and the number of correlated electrons, our aim is also presenting accurate CCSD(T) data for selected properties of this molecule. We will focus on first three electronic states, the ground 3 Σ -state, and low lying 3 Π and 1 Σ + excited states of GaN. We calculated bond lenghts, dissociation and excitation energies, dipole moments and dipole polarizabilities. The electronic configuration for the three states is ...(10σ) 2 (4π x ) 1 (4π y ) 1 , ...(10σ) 1 (4π x ) 2 (4π y ) 1 , and ...(10σ) 0 (4π x ) 2 (4π y ) 2 , respectively. GaN is an important molecule since it is a part of common semiconductor materials widely used in nanotechnology and optoelectronics. Compounds containing GaN received a lot of attention in many experimental and theoretical studies in which molecular clusters [5,[9][10][11][12] and simulations of GaN containing surface processes [13,14] were considered. Unfortunatelly, there is no experimental study of GaN diatomic at the time. A few theoretical papers

can be found in the literature. The ground state was determined as 3 Σ -. A qualitative bonding picture of this state was analyzed by Mochizuki and Tanaka [3]. Subsequent studies were presented by Ueno et al. [15], Denis and Balasubramanian [4] and recently Tzeli and Tsekouras [2]. We used a series of aug-cc-pVXZ (X=T,Q,5) and aug-cc-pCVXZ (X=T, Q, and 5) basis sets correlating three up to 21 electrons of Ga and five up to all seven electrons of the nitrogen atom. In some calculations we have correlated all 31 electrons of Ga, but actually there is no contracted basis set constructed having in mind correlating also inner shell 1s 2 2s 2 2p 6 electrons of Ga. Therefore we do not present calculations with all 38 electrons correlated for GaN in this paper. By systematic comparison of results with different basis sets and different correlation schemes we tried to reveal inaccuracies when combining a specific basis set with inappropriate number of correlated electrons. All basis sets were constructed having in mind their use within the spin-free Douglas-Kroll-Hess (DKH) relativistic approximation [8].

Along with results using the aug-cc-pVXZ (X=T,Q,5) and aug-cc-pCVXZ basis sets we also present CCSD(T)-DKH calculations with the relativistically constructed ANO-RCC (Atomic Natural Orbitals) basis set [16]. ANO-RCC basis sets for Ga and N were contracted to the size corresponding to the aug-cc-pVQZ basis. The difference is that the aug-cc-pVQZ basis basis was constructed with correlating only three valence electrons of Ga, whereas the ANO-RCC basis set was developed having in mind correlation of 13 outer shells 3d 10 4s 2 4p 1 electrons of the Ga atom.

Computational details

All calculations were carried out with the MOLCAS suite of programs [17]. The CCSD(T) Coupled Cluster method with single and double excitations and perturbative triples [18,19] with an approximate spin adaptation [START_REF] Urban | Recent Advances in Coupled-Cluster Methods[END_REF] was employed in all calculations using the Restricted Hartree Fock reference for high spin open-shell states. Our method is similar to the approach developed independently by Knowles, Hampel and Werner [START_REF] Knowles | [END_REF]. This method leads to energies which are very close to the energies obtained by the rigorous (even if slightly more demanding) method [22] for the spin adaptation. Denominators in perturbatively calculated triples were diagonal Fock-matrix elements [23], see also Refs. [START_REF] Urban | Recent Advances in Coupled-Cluster Methods[END_REF]24]. In all calculations we have considered scalar relativistic effects using the second-order spin-free scalar relativistic approach introduced by Douglas, Kroll, and Hess (DKH) [8].

The evaluation of the dipole moment and the parallel zz component of the dipole polarizability was accomplished using the numerical finite field perturbation approach [25]. The strength of the external electric field was 0.001, 0.002, and 0.004 au, respectively. Final dipole moments and dipole polarizabilities were obtained using the numerical derivatives prescription, which allow to remove the higher order terms [26,27].

Five correlation schemes were used for aug-cc-pVXZ, aug-cc-pCVXZ (X=T,Q,5), and ANO-RCC basis sets. For the gallium atom we correlated 3 valence (4s 2 4p 1 ) electrons, 13 (3d 10 4s 2 4p 1 ), 19 (3p 6 3d 10 4s 2 4p 1 ) and 21 (3s 2 3p 6 3d 10 4s 2 4p 1 ) electrons. These four correlation schemes for Ga were combined with correlating five valence electrons of the nitrogen leading to the total number of 8, 18, 24, or 26 correlated electrons of the three states of GaN. In addition, the (3s 2 3p 6 3d 10 4s 2 4p 1 ) correlation scheme of 21 electrons of Ga was combined with correlating all seven electrons of nitrogen, leading to the total number of 28 correlated electrons. Some data, not presented here, were also obtained by correlating all 31 electrons of Ga.

Basis sets used for calculations of molecular properties of GaN molecule differ in size significantly. This is important for considering the relation between the accuracy and computer demands. One should realize that the worst step in CCSD(T) scales with the number of occupied and virtual orbitals with the factor No 3 Nv 4 . Aug-cc-pVXZ basis has 105 contracted basis functions for X=T, 173 for X=Q and 267 for X=5, respectively. ANO-RCC was contracted to the size of aug-cc-pVQZ having 173 basis functions. Core-valence basis sets are considerably larger, having 159 functions in the aug-cc-pCVTZ basis set, 279 and 444 basis sets in the augcc-pCVQZ and aug-cc-pCV5Z basis sets, respectively. We note that the standard aug-cc-pV5Z and core-valence aug-cc-pCVQZ basis sets contain approximately the same number of basis functions, but were constructed differently.

In order to predict molecular properties more accurately we have estimated the complete basis set limit (CBS) using the 1/X 3 fit with respect to the basis set cardinal number X [28]. Data from aug-cc-p(C)VXZ series with X=T,Q,5 or X=Q,5 were employed. Dissociation energies were corrected for the basis set superposition error (BSSE) by the usual counterpoise correction [7].

Results

Considering the manner adopted for the construction of the aug-cc-pVXZ series, proper number of correlated electrons for GaN with these basis sets would be 8. Proper number of correlated electrons with the aug-cc-pCVXZ series is 28. Results with the aug-cc-pCVXZ series and 28 correlated electrons will be considered as reference numbers for each property. With the ANO-RCC basis set the proper number of correlated electrons is 18. Results with these "proper" contraction schemes are printed in bold in Tables 1 to 4, in which we summarize optimized CCSD(T) bond distances, BSSE corrected dissociation energies and electric properties of the GaN molecule. A transparent view upon the relation between the basis set and the number of correlated electrons is seen in Figure 1 for the excitation energy from the 3 Σ -ground state to the 3 Π state and in Figure 2 for the excitation energy to the 1 Σ + state. In Table 5 we present our best results for dissociation and excitation energies, dipole moments and dipole polariabilities for the three states of GaN with 28 correlated electrons and CBS limits obtained by extrapolations from aug-cc-pCVXZ (X=T, Q, 5) basis sets.

GEOMETRIES

Optimized bond lenghts with the aug-cc-pVXZ and ANO-RCC basis sets are presented in Table 1 and results with the aug-cc-pCVXZ series are presented in Table 2. It is immediatelly clear that results with only 8 correlated valence electrons of GaN in all three states differ significantly, by 0.10 -0.15 from R e with 28 correlated electrons and the CBS limit. This holds for the aug-cc-pVXZ, aug-cc-pCVXZ series and also for ANO-RCC basis sets reconfirming the necessity of correlating at least 3d electrons of the Ga atom containing species. To correlate 3p electrons of Ga is needed for more accurate results, remaining core-valence correlation effect is unimportant particularly for excited states. Going deeper beyond the 3d Ga shell with, say, the aug-cc-pCVQZ basis set, optimized values of R e for the ground and the two excited states of GaN are lying within the relatively narrow interval of 0.002 -0.006 bohr with 18 up to 26 correlated electrons. With the valence-only aug-cc-pVQZ basis, this interval is much larger, 0.015-0.025 Bohr leading to exagerrated importance of the core-valence correlation effects. With all aug-cc-pVXZ bases R e decreases upon correlating more than 3d 10 4s 2 4p 1 electrons of Ga, in con- trast to the trend with the aug-cc-pCVXZ series. Since aug-cc-pVXZ basis sets are not constructed for treating the core-valence correlation, increased computational effort leads in fact to inferior results. CBS limits from aug-cc-pCVXZ basis sets with 18 correlated electrons (3d 10 4s 2 4p 1 electrons of Ga) do not change significantly upon considering also the Ga 3s 2 3p 6 shell. Dominating core-valence correlation arises from the 3d correlation effect. The nitrogen core correlation effect can be estimated by comparing results with 26 and 28 correlated electrons. It amounts 0.003 -0.002 Bohr for the ground state and both excited states. The convergence of R e toward the CBS limit with the aug-cc-pCVXZ series is smooth. However, with the smallest aug-cc-pCVTZ is R e overestimated by about 0.03 bohr (more for the ground state than for excited states). Accidentally, the (Q,5) CBS limit as obtained from augcc-pVQZ and aug-cc-pV5Z basis sets agrees relatively well with the aug-cc-pCVXZ CBS limit with 24-28 correlated electrons. This is an example of a good result for a wrong reason. The basis set dependence of R e on the principal number X (X=T, Q, and 5) is oscillatory for aug-cc-pVXZ bases, with no control of accuracy, in contrast to the aug-cc-pCVXZ series, which is smooth. Should aug-cc-pVXZ basis sets be used, it is imperative to avoid the smallest member of the series, the augcc-pVTZ basis, since the (T,Q,5) extrapolated value differs from the (Q,5) CBS limit by up to 0.04 bohr. With aug-cc-pCVXZ basis sets the three-point and the two-point CBS limits agree excellently.

When thinking about the most economic way of obtaining acceptably accurate R e with a single basis set, one should be satisfied with correlating just 18 electrons in CCSD(T). For all three states of GaN good agreement with our best CBS limit is obtained with the ANO-RCC basis. It clearly outperforms the aug-cc-pVQZ basis set which has smaller number of primitives but the same number of contracted basis functions. Due to the oscillatory behaviour of R e calculated with aug-cc-pVXZ basis sets we do not recommend using aug-cc-pVTZ basis for Ga containing compounds, even if, paradoxically, result with 28 correlated electrons (3.748 bohr, as compared with the CBS limit, 3.750 bohr) is not bad, clearly thanks to the error cancellations. However, since computer demands of CCSD(T) calculations increase with the number of correlated electrons with the power of 4 the advantage of using the small basis set is diminished. Relatively small is the aug-cc-pCVTZ basis set, but due to missing cancellation of errors results are no better than with ANO-RCC or aug-cc-pVQZ basis sets. Aug-cc-pCVQZ basis set is larger but results are more reliable.

DISSOCIATION ENERGIES

The performance of the valence-only aug-cc-pVQZ, aug-cc-pCVQZ and ANO-RCC basis sets in calculations of electronic adiabatic dissociation energies (D e ) is presented in Table 3. Geometries were taken from Tables 1 and2. The ground state X 3 Σ -and the first excited state A 3 Π dissociate to Ga 2 P + N 4 S, the second excited a 1 Σ + state dissociates via the Ga 2 P + N 2 D channel.

Analyzing data in Table 3 we observe, first, that correlating just eight valence electrons provides D e s which are underestimated by about 0.06, 0.1, and 0.23 eV for the ground state and A 3 Π and a 1 Σ + excited states in comparison with calculations with 18 to 28 correlated electrons. Second, upon correlating 18 to 28 electrons all values lie within relatively narrow interval of about 0.02 -0.03 eV. This is a small value considering other effects which affect the accuracy of calculated D e . In light of the discussion on proper use of basis sets and the number of correlated electrons for calculations of the GaN bond length, dissociation energies of all three states agree very well, to within 0.02 eV or even better for aug-cc-pVQZ, aug-cc-pCVQZ, and ANO-RCC basis sets. Thus, somewhat surprisingly, there are no remarkable differences in the performance of the three basis sets in CCSD(T) calculations of D e . This holds for all three states considered in this work. Yet, there is an indication signalizing inferior performance of the aug-cc-pVQZ basis set. Good results with this basis are accompanied by quite large BSSE (see last three columns of Table 3), 0.1 up to 0.4 eV for the three states. The only exception is D e calculated with 8 correlated electrons. Large BSSE compensates the defficiency of the insufficient description of the deeper electron shells of the Ga atom with the aug-cc-pVQZ basis set. BSSE is small for ano-rcc basis when we correlate 3d electrons of gallium (which is exactly what it was constructed for). In case of aug-cc-pCVQZ basis BSSE is order of magnitude smaller than with the aug-cc-pVQZ basis for all correlation schemes. Exaggerated importance of the BSSE correction when the basis set is not appropriate for the number of electrons which are correlated was also observed in calculations of D e for the ground state of GaN and for Ga 2 by Tzeli and Tsekouras [2]. Their CCSD(T) CBS limit for D e of the ground state (2.12 eV) obtained with 18 correlated electrons employing small-core relativistic pseudopotentials basis agrees reasonably well with our best value, 2.140 eV (Table 5).

As mentioned above, the second excited a 1 Σ + state dissociates via the Ga 2 P + N 2 D channel. Since the doublet state of the nitrogen atom shows a multireference character which is accompanied by high excitation amplitudes (up to 0.28) in the single reference CCSD calculation, results with the single-determinant reference CCSD(T) method may not be conclusive due to quasidegeneracy effects. To ensure that our conclusions concerning the performance of differently constructed basis sets are not plagued by quasidegeneracy effects, we used the CASPT2 method for evaluation of dissociation energies of all three electronic states along with CCSD(T). The active space contains 8 valence electrons in 8 active orbitals (4s4p of the gallium atom and 2s2p orbitals of the nitrogen atom). CASPT2 dissociation energies with 26 correlated electrons using aug-cc-pVQZ, aug-cc-pCVQZ and ANO-RCC basis sets are presented in Table 3. Comparison of CCSD(T) and CASPT2 dissociation energies with the series of aug-cc-pCVXZ basis sets (X=T, Q. and 5) is shown in Table 5 for 28 correlated electrons. For the ground X 3 Σ -state CCSD(T) and CASPT2 results differ by less than 0.1 eV, CASPT2 D e s values being sligtly higher than CCSD(T) dissociation energies. Similarly, CASPT2 dissociation energies of the first A 3 Π excited state are higher than CCSD(T) ones, this time by 0.12 to 0.17 eV. Conversely, for the a 1 Σ + state is the CASPT2 dissociation energy lower than the CCSD(T) energy. This reflects the overestimated energy of the nitrogen 2 D state by CCSD(T). The importance of triples in CCSD(T) can be determined by comparing CCSD(T) and CCSD dissociation energies in Table 5. For the ground state this contribution is relatively small, 0.1 eV, but for the A 3 Π and the a 1 Σ + states it increases up to 0.2 and almost 0.5 eV, respectively. Therefore, comparisons of CCSD(T), CCSD, and CASPT2 data demonstrate that details of the treatment of correlation effects matter. Concerning BSSE, conclusions from CASPT2 and CCSD(T) data for D e are completely the same: BSSE is much larger for aug-cc-pVQZ basis set than for the aug-cc-pCVQZ basis. With CASPT2 it is larger than with CCSD(T) with all three bases.

DIPOLE MOMENTS AND DIPOLE POLARIZABILITIES

The performance of aug-cc-pVQZ, aug-cc-pCVQZ and ANO-RCC basis sets in calculations of dipole moments and parallel dipole polarizabilites with 8 to 28 correlated electrons is documented in

are taken from Tables 1 and2 for each particular correlation scheme and the basis set. Originally, we expected that dipole moments and polarizabilities will be very sensitive to the number of correlated electrons and the basis set. It appears, however, that the sensitivity of electric properties to the basis set is not dramatic. For example, considering 26 correlated electrons for the ground state with the augcc-pVQZ basis set (inappropriate basis for this calculation) the dipole moment is 0.719 a.u. It differs by 9.1 % from the aug-cc-pCVQZ basis set result. Note that ANO-RCC result differs from this reference number much less, by 3.3 %. The basis set dependence for both dipole moment and polarizability with the aug-cc-pCVXZ series and 28 correlated electrons is even smaller, as expected, see Table 5.

Aside from data obtained with 8 correlated electrons, dipole moments and dipole polarizabilities are insensitive also to the number of correlated electrons. With all basis sets the CCSD(T) dipole moments of GaN with 18-26 correlated electrons lie in the interval narrower than 0.02 a.u. for all states. Core correlation effect due to correlating inner shell 1s electrons of nitrogen is 0.005 a.u. for the a 1 Σ + state and even smaller for the ground and the A 3 Π excited states. Analogous insensibility holds for the parallel polarizability.

EXCITATION ENERGIES

Adiabatic CCSD(T)-DKH excitation energies T e from the ground X 3 Σ -state to the first two excited A 3 Π and a 1 Σ + states, respectively, are much more sensitive to the selection of the basis set and proper correlation scheme than are other properties. Nevertheless, when disscussing details of the basis set effects we should realize that the energy gap between the ground X 3 Σ -state and the A 3 Π excited state is quite small. For the two states with 8-28 correlated electrons the basis set dependence of T e is transparently demonstrated in Figures 1 and2, respectively. On the left hand side we present the series aug-cc-pVXZ (X= T, Q, and 5), supplemented with ANO-RCC results. On the right hand side there are excitation energies with the aug-cc-pCVXZ series. For both series we also present CBS limits. Not much discussion is needed for results with 8 correlated electrons. Excitation energies with all basis sets are quite similar, the basis set dependence for both aug-cc-pVXZ and aug-cc-pCVXZ series is smooth, but excitation energies are much too high. More interesting are results with 18-28 electrons. As a curiosity we note that T e with the smallest aug-cc-pVTZ basis and 18 correlated electrons agrees very well with the CBS limit calculated from much more demanding aug-cc-pCVXZ series, 0.104 and 0.747 eV for A 3 Π and a 1 Σ + states, respectively. Typically, this is a good result for the wrong reason. Upon correlating 26 or 28 electrons and extending the aug-cc-pVTZ basis to aug-cc-pVQZ the excitation energy of the A 3 Π state diminishes to about 0.04 eV. Extrapolating to the CBS limit with these two basis sets we would obtain completely misleading result, 0.007 eV with 28 correlated electrons. This may lead to some doubts about the state ordering in GaN. We remind that there was a long lasting discussion about the ordering of states for valence isoelectronic BN and AlN molecules stressig the importance of proper handling of electron correlation effects for these species [29,30]. When calculating the CBS limit using smaller aug-cc-pVXZ series one should be careful since the behavior of T e with respect to the cardinal number X is oscillatory for both excited states. In any case, the smallest aug-cc-pVTZ must be avoided in the extrapolation to the CBS limit when core electrons are correlated. In contrast, T e is much less sensitive to the cardinal number in the aug-cc-pCVXZ series, it is also almost insensitive to the number of correlated electrons, the basis set dependence is smooth and thus results are much more reliable. The problem is the size of aug-cc-pCVXZ bases. The largest aug-cc-pCV5Z basis set with 444 basis functions is no problem for the CCSD(T) calculation of a single GaN molecule but would be demanding or intractable for dimers or larger complexes. For the sole calculation of T e a good choice is using the ANO-RCC basis set. Another good candidate for a single-shot calculation would be the aug-cc-pCVQZ basis set (279 basis functions). For the sake of saving the computer time with this basis set it is satisfactory to correlate only 18 electrons (even if it was constucted having in mind 28 correlated electrons for GaN). Much larger error for both excited states in comparison with the CBS limit from Table 5, is obtained with, e.g., the aug-cc-pV5Z basis set (267 basis functions).

For accurate calculation of excitation energies of GaN and other IIIA element molecules it is necessary to consider carefuly also the treatment of the electron correlation. Starting from InN one should also take care of the spin-orbit coupling [31]. We observe (Table 5) that T e for the A 3 Π state with the CCSD method is higher than with CCSD(T) by about 0.1 eV and even more (by as much as 0.5 eV) for the a 1 Σ + state. Therefore, triples contribute significantly. At the same time, with CASPT2 (8,8) are excitation energies for the A 3 Π state lower than CCSD(T) energies by about 0.06-0.08 eV with aug-cc-pCVXZ basis sets. Finally, users of the valence-only aug-cc-pVXZ basis sets should have been alert to the possibility of the reversed order of stabilities of states with CASPT2. With the aug-cc-pVQZ basis set T e of the A 3 Π state decreases with the number of correlated electrons from 0.154 eV for 8 electrons correlated even up to the negative value of -0.004 eV for 26 electrons correlated.

CONCLUSIONS

Our best CCSD(T) DKH results for dissociation and excitation energies, dipole moments and dipole polarizabilities of the 3 Σ -, 3 Π and 1 Σ + states of the GaN molecule are summarized in Table 5. All results include the scalar relativistic effects and the core-valence electron correlation effects treated with the aug-cc-pCVXZ basis set series. Our CBS limits can be considered as benchmark results for properties of the low-lying states of GaN. The CCSD(T) DKH dissociation energy of all three states of GaN should be supplemented by the spin-orbit contribution. According to recent work of Shirai [33] the value for the 2 P 3/2 -2 P 1/2 splitting of the Ga atom is 829 cm -1 . The spin-orbit splitting for the 3 Σ -state of the GaN molecule is negligible [4,32]. Final best estimate of D e for the ground 3 Σ -state is 2.072 eV. Our CCSD(T) DKH dissociation energy as presented in Table 5 agrees reasonably well with value (2.016 eV) published by Denis and Balasubramanian [4]. Along with CCSD(T) they employed also the MRSDCI method. Similarly as it follows from our CASPT2 and CCSD(T) results, MRSDCI method predicts a smaller separation between the 3 Σ -and the 3 Π states, 532 cm -1 in comparison with CCSD(T) (857 cm -1 ). Even if they do not provide the effect of the basis set in relation to the core-valence correlation their results are in qualitative agreement with our conclusions. Both CCSD(T) and MRSDCI methodologies indicate that the ground state of GaN is 3 Σ -. Our results also compare very well with more recent results by Tzeli and Tsekouras [2]. These authors studied several gallium and nitrogen containing molecules including the X 3 Σ -state of GaN. They correlated 8 and 18 electrons, respectively, using small-core pseudopotentials and various basis sets. The difference between our best geometries and their results for the ground state is only 0.001 Bohr, dissociation energies differ by only 0.02 eV. General conclusions and trends concerning the performance of aug-cc-pVXZ and aug-cc-pCVXZ basis sets (with X going up to quintuple zeta for both series in the present study) in accurate calculations of these properties of the GaN ground X 3 Σ -state are the same. These result are supplemented in this paper by correlating deeper shells of GaN and its low-lying excited states. It appears that excitation energies are most sensitive to proper selection of the basis set as used with different correlation schemes. Our excitation energy for the A 3 Π state agrees reasonably well with the value published by Denis and Balasubramanian [4], 0.106 eV. However, their value for the a 1 Σ + state is overestimated in comparison with our CBS limit (Table 5) by 0.05 eV. This state is more sensitive to basis set effects. Concerning electric properties, we were only able to compare electric dipole moments from the work of Denis and Balasubramanian, obtained by the multireference MRSDCI method finding the differences from our results of about 0.1 a.u. for the three states.

Our systematic investigation of the accuracy, reliability and efficiency for calculations of molecular properties of GaN as a typical species containing a IIIA element leads us to several conclusions on the use of differently constructed basis sets for the treatment of the core-correlation effects. Clearly, it would be advantageous to use relatively smaller aug-cc-pVXZ basis sets for considering the core-valence correlation effects. These basis sets are usually flexible enough in many applications and are smaller than their aug-cc-pCVXZ counterparts. Clearly, some data obtained with these basis sets agree with the CBS limit quite well. However, we were unable to obtain good results in a systematic and controlled manner. Simply, these basis sets were not constructed for core-valence correlations. Yet, considering at least 3d electrons for such species is inevitable. In light of this fact and the unreliable performance of the valence aug-cc-pVXZ basis sets when correlating deeper shels of Ga we see little space for using these basis sets as a predictive tool for properties of molecules containing these atoms. Main obstacle is our finding that molecular properties calculated with aug-cc-pVXZ basis sets suffer from the oscillatory behavior with increasing cardinal number. Smaller members of the series, as the aug-cc-pVTZ basis set, should be avoided in any extrapolations to the CBS limit when core electrons are correlated. Then the advantage of its relatively smaller size is lost. For larger-scale calculations of IIIA elements compounds we would prefer smaller members of the aug-cc-pCVXZ series. A typical example is the aug-cc-pCVQZ basis set (279 basis functions) which mostly outperforms the aug-cc-pV5Z basis with 267 basis functions. Large number of virtuals can be circumvented by using the truncated space obtained by the OVOS (Optimized Virtual Orbital Space) technique [34]. Utilizing the fact that for virtually all properties calculated in this paper the effect of core valence electrons beyond 3d was relatively small we can restrict ourselves just to these electrons. For a single basis set calculation very good compromise between the accuracy and efficiency offer the ANO-RCC basis sets.

F

. CCSD(T)-DKH BSSE corrected dissociation energies (eV) and BSSE (eV) of 3 Σ -b , 3 Π b and 1 Σ + c electronic states of GaN employing various basis sets. Geometries of GaN are taken from Tables 1 and2. Data in parenthesis with 26 correlated electrons were calculated using the CASPT2 (8,8)

Molecular Physics GaN˙rev

Table 5.

Excitation energies (eV), BSSE corrected dissociation energies a (eV), dipole moments (a.u.) and polarizabilites (a.u.) of first tree states of the GaN molecule obtained at CCSD(T) DKH level using aug-cc-pCVXZ (X=T,Q,5) basis sets and correlating 28 electrons. Two extrapolation schemes employ the three-point aug-cc-pCVXZ (X=T,Q,5) data (the TQ5 extrapolation) and the two-point aug-cc-pCVXZ (X=Q,5) data (the Q5 extrapolation). Geometries used are the (Q,5) extrapolated values from Table 2. Numbers in parentheses are CASPT2 (8,8)

Figure 2 .

 2 Figure 2. CCSD(T) DKH excitation energies for GaN a 1 Σ + state. The left hand side: Series of valence aug-cc-pVXZ basis sets, with empty symbols represening ANO-RCC basis set contracted to the size of the aug-cc-pVQZ basis. Right hand side: series of core-valence aug-cc-pCVXZ basis sets.

Table 4 .

 4 Optimized internuclear distances

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1 .

 1 Bond lenghts (Bohr) of the 3 Σ -, 3 Π and 1 Σ + electronic state of GaN molecule obtained at the CCSD(T) DKH level employing series of valence basis sets with ANO-RCC basis set and usig 5 correlation schemes.

	o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 2 .

 2 Bond lenghts (Bohr) of the 3 Σ -, 3 Π and 1 Σ + electronic state of GaN molecule obtained at CCSD(T) DKH level employing series of core -valence basis sets and usig 5 correlation schemes.

		No. of		Basis set and its size		Extrapolation
	State	corr. el.	aug-cc-pCVTZ	aug-cc-pCVQZ	aug-cc-pCV5Z	(T,Q,5)	(Q,5)
			159	279	444
	X -3 Σ -	8	3.866	3.858	3.857	3.854	3.856
		18	3.781	3.762	3.755	3.748	3.748
		24	3.788	3.768	3.761	3.754	3.754
		26	3.787	3.767	3.760	3.753	3.753
		28	3.784	3.764	3.757	3.750	3.750
	A -3 Π	8	3.621	3.614	3.613	3.610	3.612
		18	3.508	3.491	3.484	3.478	3.477
		24	3.515	3.497	3.489	3.482	3.481
		26	3.514	3.496	3.489	3.482	3.482
		28	3.512	3.494	3.487	3.480	3.480
	a -1 Σ +	8	3.328	3.323	3.322	3.320	3.321
		18	3.221	3.209	3.204	3.200	3.199
		24	3.223	3.210	3.205	3.200	3.200
		26	3.224	3.211	3.206	3.201	3.201
		28	3.222	3.209	3.204	3.199	3.199

Table 4 .

 4 method. -basis set was contracted to size aug-cc-pVQZ basis set. b -dissociation products were 2 P Ga and 4 S N. c -dissociation products were 2 P Ga and 2 D N. CCSD(T)-DKH dipole moments (a.u.) and parallel dipole polarizabilites (a.u.) of the 3 Σ -, 3 Π and 1 Σ + electronic states of the GaN molecule with aug-cc-pVQZ, aug-cc-pCVQZ and ANO-RCC basis sets and 8-28 correlated electrons. Geometries taken from Tables1 and 2. -basis set was contracted to size aug-cc-pVQZ basis set.

		No. corr.	Dissociation energies (eV)		BSSE (eV)
		electrons	aug-cc-pVQZ	aug-cc-pCVQZ	ano-rcc a	aug-cc-pVQZ	aug-cc-pCVQZ	ano-rcc a
	X -3 Σ -	8	1.999	2.009	2.002	0.007	0.005	0.007
		18	2.062	2.071	2.065	0.133	0.015	0.045
		24	2.048	2.061	2.055	0.198	0.018	0.104
		26	2.050 (2.156)	2.064 (2.146)	2.058 (2.147)	0.220 (0.258)	0.018 (0.037)	0.132 (0.145)
		28	2.051	2.065	2.059	0.230	0.019	0.217
	A -3 Π	8	1.824	1.832	1.827	0.008	0.006	0.008
		18	1.949	1.959	1.956	0.173	0.018	0.051
		24	1.933	1.946	1.944	0.264	0.022	0.118
		26	1.934 (2.058)	1.949 (2.105)	1.947 (2.113)	0.294 (0.359)	0.023 (0.042)	0.150 (0.162)
		28	1.935	1.949	1.948	0.305	0.023	0.244
	a -1 Σ +	8	3.733	3.750	3.751	0.010	0.008	0.010
		18	3.960	3.966	3.972	0.221	0.024	0.060
		24	3.925	3.932	3.939	0.345	0.028	0.137
		26	3.934 (3.749)	3.944 (3.788)	3.951 (3.811)	0.384 (0.461)	0.029 (0.051)	0.173 (0.182)
		28	3.955	3.970	3.974	0.397	0.030	0.276

a a

  results, see text. -spin-orbit correction for Ga 2 P is 0.068 eV. b -used field strength 0.002 0.004 and 0.008.

				Basis set		Extrapolation
	State	Property	aug-cc-pCVTZ	aug-cc-pCVQZ	aug-cc-pCV5Z	(T,Q,5)	(Q,5)
	X -3 Σ -	De	1.960(2.023)	2.064(2.147)	2.101(2.199)	2.140(2.244)	2.140(2.254)
		De(CCSD)	1.836	1.933	1.966
		µz	0.773	0.780	0.781	0.784	0.782
		αzz	54.23	53.63	53.42	53.20	53.20
		αxx	50.72 b	49.74 b
	A -3 Π	Te	0.122(0.063)	0.112(0.038)	0.108(0.030)	0.104	0.104
		Te(CCSD)	0.221	0.214	0.212
		De	1.827(1.934)	1.949(2.103)	1.991(2.155)	2.037(2.219)	2.035(2.210)
		De(CCSD)	1.606	1.715	1.753
		µz	1.065	1.074 b	1.076 b	1.080	1.078
	a -1 Σ +	αzz Te Te(CCSD) De De(CCSD) µz αzz F o 0.829(0.752) 57.48 1.308 3.833(3.550) 3.391 1.542 54.16 r F o r	57.24 b 0.780(0.702) 1.274 3.969(3.793) 3.502 1.547 53.81	57.06 b 0.764(0.686) 1.262 4.015(3.835) 3.539 1.548 53.63	56.98 0.746 4.066(3.931) 1.550 53.51	56.87 0.747 4.063(3.879) 1.549 53.44
			P P
			e e
			e r e r
				R R
				e e
				v i e v i e
					w w
					O n l O n l
						y y

a

Michal Kraus, Lucia Šimová, Pavel Neogrády and Miroslav Urban

Acknowledgements

This work was supported by the Slovak Research and Development Agency, contracts No. APVV-20-018405, APVV-LPP-0155-09, and the Slovak Grant Agency VEGA under the contract No. 1/0428/09. Support from EURATOM, contract No.FU07-CT-2006-00441 is also gratefully acknowledged.


H W Mensink
email: h.mensink@oogziekenhuis.nl

J Vaarwater

R J W De Keizer

D De Wolff-Rouendaal

C M Mooy

A De Klein

D Paridaens

Chromosomal aberrations in iris melanomas

Keywords: iris melanoma, FISH, monosomy 3, CDKN2A

Background. Uveal melanomas can develop in the choroid, ciliary body and in the iris. In choroidal and ciliary body melanomas specific chromosomal changes correlate with metastatic disease. Iris melanomas have a better prognosis than choroidal melanomas and it would be interesting to know if they share chromosomal changes. In addition, iris melanomas might harbour UV-induced mutations of tumour suppressor genes, such as PTEN and CDKN2A.

Methods. Twenty iris melanomas were analyzed for chromosome 1p, 3, 6, 8, 9p and 10q abnormalities using fluorescence in situ hybridization. These results were correlated to clinical follow-up data using statistical analyses.

Results. (Partial) loss of chromosome 3 was observed in 9 iris melanomas and gain of 8q was present in 7 tumours. Loss of chromosome 9p was demonstrated in 7 tumours but no deletions of the PTEN region on chromosome 10 were found. Three patients died of metastatic disease and one patient developed liver metastases, but is still alive. Univariate analysis indicated a lower disease free survival for patients with diffuse growing melanomas (p=0.01), melanomas that lost a copy of chromosome 3 (p=0.03), or invading the ciliary body (p=0.01). In multivariate analysis none of the correlations were significant.

Conclusion. Loss of chromosome 3 as well as loss of chromosomal region 9p21 (that entails tumour suppressor gene CDKN2A) play a role in iris melanoma. A firm correlation with disease-free survival could not be established, possibly due to the small sample size,

Introduction

Iris melanomas are the least common of uveal melanomas (< 6%) and tend to present at a smaller size than posterior melanomas. [START_REF] Henderson | Iris melanoma[END_REF] They can cause drainage blockage of the anterior chamber angle and secondary elevation of the intraocular pressure. They carry a better prognosis compared to other uveal melanomas with a metastatic rate of 0.5% at 3 and 7% at 10 years, respectively. [START_REF] Shields | Metastasis of uveal melanoma millimeter-by-millimeter in 8033 consecutive eyes[END_REF] Melanocytic lesions of the iris are usually managed by a 'wait and see policy' with regular follow-ups to detect possible growth. Biomicroscopy, gonioscopy and measurements with the slit lamp are standard tests. Ancillary tests such as fluorescein angiography and ultrasound biomicroscopy are sometimes needed to assess growth more accurately. Iris melanomas can be treated by local excision, enucleation or radiotherapy. Metastases are more likely to develop in those patients who are older or with diffuse iris melanomas, compared to well-defined nodular ones. Other warning signs are tumour features of iris root location and ciliary body invasion with elevated intraocular pressure and extraocular extension. [START_REF] Shields | Iris melanoma: risk factors for metastasis in 169 consecutive patients[END_REF] The differential diagnosis of iris melanoma includes Cogan Reese syndrome, nevus, melanocytoma, metastasis, adenoma, adenocarcinoma, cysts and pigment dispersion syndrome. [START_REF] Shields | The differential diagnosis of malignant melanoma of the iris. A clinical study of 200 patients[END_REF] The favourable outcome of iris melanomas compared to choroidal and ciliary body melanomas is poorly understood. To determine the prognosis on only histological classification proved to be difficult. [START_REF] Jakobiec | Are most iris "melanomas' really nevi? A clinicopathologic study of 189 lesions[END_REF] Therefore, more accurate tests are needed. The different anatomical locations of the intraocular melanomas result in different microenvironments. Aqueous humour and vitreous contain various soluble factors that might influence cell growth. [START_REF] Mudhar | The in vivo modulatory effects of an anterior-chamber microenvironment on uveal melanoma[END_REF] Another differential factor is UV-radiation. Its role in uveal melanoma is inconclusive, but arguably the iris is more exposed to the sun than the choroid. [START_REF] Singh | Sunlight exposure and pathogenesis of uveal melanoma[END_REF][START_REF] Mooy | No N-ras mutations in human uveal melanoma: the role of ultraviolet light revisited[END_REF][START_REF] Schmidt-Pokrzywniak | Positive interaction between light iris color and ultraviolet radiation in relation to the risk of uveal melanoma: a case-control study[END_REF][START_REF] Li | Patterns of tumor initiation in choroidal melanoma[END_REF][START_REF] Manning | Ocular melanoma in a Long Evans rat[END_REF][START_REF] Marshall | The effect of blue light exposure and use of intraocular lenses on human uveal melanoma cell lines[END_REF][START_REF] Vajdic | Sun exposure predicts risk of ocular melanoma in Australia[END_REF] While UV-radiation might induce changes in the genome, intrinsic genetic differences may also play a role. [START_REF] Sisley | Cytogenetics of iris melanomas: disparity with other uveal tract melanomas[END_REF] To further elucidate the genetic differences between choroidal and iris melanomas, the following study was performed. Twenty iris melanomas were analyzed for chromosomal aberrations well known in choroidal and ciliary body melanomas, such as loss of chromosome 3 and gain of chromosome 8q [START_REF] Singh | Uveal melanoma: genetic aspects[END_REF] , in addition to common changes in skin melanomas, such as loss of 9p21 and 10q23 (regions that harbour the tumour suppressor genes CDKN2A and PTEN, respectively). [START_REF] Wu | PTEN signaling pathways in melanoma[END_REF]

Materials and methods

The research methodology followed the tenets of the Declaration of Helsinki. Informed consent was obtained for the use of human tissue from the patients. The research was approved by the local research ethics committee. All cases of iris melanomas were collected from the files from the departments of pathology in Rotterdam (1993-2009) and Leiden (1980 to 2009) and reviewed retrospectively. Cases that originated from the ciliary body were excluded as well as patients who had no available follow-up information. This resulted in 22 iris melanomas available for our study. FISH (Fluorescence in situ Hybridisation) was unsuccessful in 2 cases.

The enucleation and biopsy specimens were fixed in standard 10%, buffered formalin for 24 h. They were embedded in conventional, paraffin wax and sections were cut at 5 µm.

Sections were stained with Haematoxylin and Eosin (H&E) in order to obtain a histopathological diagnosis. The specimens were reviewed and classified according to the Jakobiec histopathologic classification of iris melanocytic tumours 5 : group 1 melanocytosis, group 2 melanocytoma, group 3 epithelioid cell nevus, group 4 intrastromal spindle cell nevus, group 5 spindle cell nevus with surface plaque, group 6 borderline spindle cell nevus, group 7 spindle cell melanoma, group 8 mixed spindle and epithelioid cell melanoma, group 9 epithelioid cell melanoma.

Paraffin-embedded iris melanoma slides were used for FISH (chromosome 1p, 3, 6, 8, 9p and 10q). The pathologist located the tumour and dual colour FISH (digoxygenin and biotin labelled probes) was performed with the following probes: RP11-48E9 (1p36), RP11-522N9 (3p13), Pα3.5 (centromere 3), RP11-64F6 (3q25), RP11-356B3 (6p22), RP11-787I22 (6q21), RP11-24P4 (8p21), D8Z2 (centromere 8), RP11-88J22 (8q22), P1.063 (9p21) and PAC190D6 (10q23). [START_REF] Naus | Detection of genetic prognostic markers in uveal melanoma biopsies using fluorescence in situ hybridization[END_REF] The concentration for centromeric probes was 5 ng per slide and for BAC probes 50 to 75 ng per slide was used. As a control for aneuploidy a probe on chromosome 5 was used: RP11-1059N10 (5q12). The tissue sections were pretreated by dewaxing with xylene, permeabilizing with sodium thyocyanate, proteolysis and denaturation. After hybridization and washing, slides were counterstained with 4', 6-diamidino-2-phenylindole and mounted in anti-fade solution (Dabco-Vectashield 1:1). The amount of scored nuclei depended on the size of the tumour, at least 50 nuclei and till a maximum of 200 nuclei were scored. Cut-off limits for deletion (>25% of the nuclei with one signal) or amplification (>10% of the nuclei with 3 or more signals) were adapted from own experiments and the available literature. [START_REF] Van Dekken | Cytogenetic analysis of human solid tumors by in situ hybridization with a set of 12 chromosome-specific DNA probes[END_REF] The primary endpoint for disease free survival (DFS) was the time to development of metastatic disease, whereby death due to other causes was being censored. Statistical analyses were performed with SPSS software, release 15.0. Actuarial probabilities of disease free survival (with an event defined as development of metastatic disease or metastasis-related death) were estimated according to the Kaplan-Meier method. To examine the possibility that other clinical, histopathological or chromosomal variations affected the prognosis, we performed Cox proportional hazard analysis for each confounding variable. An effect was considered significant if the p-value was 0.05 or less. The following factors were evaluated: age, cell type (spindle, mixed or epithelioid cell), Jakobiec classification, growth pattern (nodular, multiple nodules or diffuse), involvement of ciliary body, tumour invasion of the chamber angle, chromosomal aberrations (chromosomes 1, 3, 6p, 6q, 8q, 9p and 10q).

Results

Tumour characteristics

Twenty iris melanomas were included in our study. The specimens were classified according to the Jakobiec histopathologic classification of iris melanocytic tumours 5 : 3 borderline spindle cell nevi, 4 spindle cell melanomas, 6 mixed spindle and epithelioid cell melanomas and 7 epithelioid cell melanomas (Figure 1).

Fifteen melanomas grew in the chamber angle resulting in secondary glaucoma in 5 cases caused by angle closure. Nine iris melanomas that spread into the angle also spread into the ciliary body, but primarily originated from the iris and in all of them epithelioid cells were present (Table 1).

Clinical data

The median age of the patients at the time of diagnosis was 62 years. Three patients had a trabeculectomy for glaucoma in their medical history (melanoma # 2, 9 and 14). Surgery or Ruthenium plaque irradiation was chosen to treat the iris melanoma.

The median follow-up time was 42 months (range 2-204 months). Three patients died of metastases. One patient had been treated with a trabeculectomy 12 years earlier (histopathological diagnosis at that time: Cogan Reese syndrome) and died of metastatic disease (melanoma # 9). The second patient presented with metastases, located in the bones and submandibular lymph nodes, 3 years after a trabeculectomy (histopathological diagnosis at that time: iris nevus without signs of malignancy) (melanoma #1). The third patient died of liver metastases 3 years after the iris melanoma, treated by Ruthenium plaque irradiation, recurred. This recurrence was treated by enucleation (melanoma #13).

Another patient presented with liver metastases 2,3 years after enucleation as the primary treatment of her iris melanoma (melanoma #10).

FISH results

All 20 iris melanoma tissue slides could be hybridized with the chromosome 3p and centromere 3 probe. For the other FISH probes it was not always possible to obtain data, because often the material was too sparse to give a definite answer. See table 2 for the results.

FISH Analysis of paraffin sections demonstrated (partial) loss of chromosome 3 in nine tumours. Seven melanomas with loss of chromosome 3 had concomitantly gained chromosome 8q (see Figure 2). Loss of the chromosomal region 9p21 was present in seven tumours, and six of those had concomitant loss of chromosome 3 and 9p. Five tumours had all three chromosome aberrations (loss of chromosome 3, 9p and gain of 8q) mentioned above. Chromosome 10q was normal in all cases.

The tumours of the patients who died of metastases revealed loss of chromosome 3. Two of them also lost 9p21, and two melanomas had gained 8q. Two of the progressively growing melanomas were composed of mixed spindle and epithelioid cells and one was composed of epithelioid cells only. All three grew into the angle of the chamber and into the ciliary body.

Three melanomas in our group were spindle cell nevi with borderline malignancy (Jakobiec group 6). Interestingly, two of them did not have any chromosomal aberrations, but the third one did lose one copy of chromosome 3q and 9p simultaneously with gain of chromosome 8q.

Statistical analysis

Univariate analysis was performed for all parameters, showing a lower disease free survival for patients with diffuse growing melanomas (p=0.01), melanomas that lost a copy of chromosome 3 (p=0.03) or invading the ciliary body (p=0.01). The other chromosomal aberrations did not show significant correlations. In a multivariate Cox proportional hazard analysis including growth pattern, ciliary body involvement and chromosome 3, they all lost significance.

Discussion

Our study reports on chromosomal abnormalities in the largest group of iris melanomas described thus far. We demonstrated that in iris melanomas (partial) loss of chromosome 3 tends to occur often, in 9/20 tumours (45%), but only correlates with progressive disease in a univariate analysis. The frequency of chromosome 3 loss is almost similar to choroidal melanomas, in which case it strongly associates with patients at risk of metastatic disease. [START_REF] Prescher | Prognostic implications of monosomy 3 in uveal melanoma[END_REF] Similar to common changes in skin melanomas, such as loss of 9p21 and 10q23, 7/ 20 iris melanomas (35%) revealed loss of chromosome 9p. UV-radiation induces mutations in the CDKN2A gene that encodes for the p16 protein that acts as a negative regulator of the cell cycle and is located on chromosome 9p21. [START_REF] Hussein | Ultraviolet radiation and skin cancer: molecular mechanisms[END_REF] Alterations of chromosome 9 occur frequently in cutaneous melanomas and relate to progressive disease. [START_REF] Wang | Constitutional alterations in p16 in patients with uveal melanoma[END_REF][START_REF] Hoglund | Dissecting karyotypic patterns in malignant melanomas: temporal clustering of losses and gains in melanoma karyotypic evolution[END_REF][START_REF] Nelson | Chromosome abnormalities in malignant melanoma: clinical significance of nonrandom chromosome abnormalities in 206 cases[END_REF] . Conversely, CDKN2A mutations are rarely observed in UM. [START_REF] Naus | Mutation analysis of the PTEN gene in uveal melanoma cell lines[END_REF][START_REF] Singh | Familial uveal melanoma: absence of germline mutations involving the cyclin-dependent kinase-4 inhibitor gene (p16)[END_REF] Our results indicate a difference between iris and posterior melanomas, because loss of 9p is found in 7/20 iris melanomas.

Yet, we could not demonstrate aberrations of chromosome 10q in iris melanomas.

Previously, it has been reported that the tumour suppressor gene PTEN (located on chromosome 10q) is not involved in uveal melanomas, whereas mutations of PTEN are present in up to 40% of cutaneous melanomas. [START_REF] Naus | Mutation analysis of the PTEN gene in uveal melanoma cell lines[END_REF][START_REF] Singh | Familial uveal melanoma: absence of germline mutations involving the cyclin-dependent kinase-4 inhibitor gene (p16)[END_REF][START_REF] Tsao | Identification of PTEN/MMAC1 alterations in uncultured melanomas and melanoma cell lines[END_REF] The same holds true for N-ras mutations. Mutations of the N-ras gene, located on chromosome 1p13, correlate with sunlight exposure and occur frequently in cutaneous melanoma arising in sun exposed areas. [START_REF] Singh | Sunlight exposure and pathogenesis of uveal melanoma[END_REF] Nevertheless, N-ras mutations are usually absent in uveal melanomas, further questioning the role of sunlight in the pathogenesis of uveal melanoma. [START_REF] Mooy | No N-ras mutations in human uveal melanoma: the role of ultraviolet light revisited[END_REF][START_REF] Singh | Uveal melanoma: genetic aspects[END_REF] In 8/20 iris melanomas the ciliary body was invaded by tumour cells, and four of these tumours spread systemically. Ciliary body involvement was only found to be a significant prognostic parameter in univariate analysis.

Little is known about somatic mutations and chromosomal aberrations in iris melanomas.

To date only 8 iris melanomas have been cytogenetically analyzed and described worldwide. The first karyogram of an iris melanoma revealed extra copies of chromosomes 2, 7 and 18, and loss of chromosome 15. Additionally, rearrangements on chromosomes 8q, 4q, 12p and marker chromosomes were observed. [START_REF] White | Cytogenetic characterization of an iris melanoma[END_REF] Sisley et al. karyotyped 3 iris melanomas and did not find specific chromosomal changes common to all tumours. Clonal abnormalities of chromosomes 3, 5, 6, 7, 8, 9, 12, 15, 17, 18, 19, and Y were found, and in one case a large number of marker chromosomes were observed. [START_REF] Sisley | Cytogenetics of iris melanomas: disparity with other uveal tract melanomas[END_REF] They concluded that iris melanomas experience relatively high levels of chromosomal alterations compared to choroidal melanomas, which is perhaps reminiscent of cutaneous melanoma. The third study included 4 iris among other uveal melanomas that were analyzed by comparative genetic hybridization. [START_REF] Vajdic | Chromosomal gains and losses in ocular melanoma detected by comparative genomic hybridization in an Australian population-based study[END_REF] Gains of chromosome 6p and 8q were identified twice and gains of chromosome 8p and 13q were identified once. One melanoma was completely diploid. These three studies did not identify critical chromosomal regions and, furthermore, the results were not compared to disease free survival.

In contrast, we observed involvement of chromosome 9p in seven tumours and loss of chromosome 3 in nine iris melanomas. These results suggest that the tumour suppressor gene CDKN2A, located on chromosome 9p, might play a role in iris melanomas, and requires analyses of the gene itself. The role of chromosome 3 loss in iris melanomas needs further investigation. The prognostic significance is not as strong as in choroidal melanomas, but possibly this is the result of the small sample size. The differences between choroidal and iris melanomas suggest that environmental influences such as UVradiation, aqueous and vitreous humour may contribute to the metastatic potential of the tumour.

The Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence on a worldwide basis to the BMJ Publishing Group Ltd and its Licensees to permit this article (if accepted) to be published in BJO editions and any other BMJPGL products to exploit all subsidiary rights, as set out in our licence (http://group.bmj.com/products/journals/instructionsfor-authors/licence-forms).

Figure 1 .

 1 Figure 1. Histology (H&E stain) of A and C) shows a spindle cell iris melanoma, and B and D) no. 1, shows a epithelioid cell iris melanoma with invasion of the angle chamber .

Figure 2 .

 2 Figure 2. FISH analysis of iris melanoma no.1, 149, demonstrates 2 signals for the biotin probe (red signal) on chromosome 10q and 1 signal for the digoxenin probe (green signal) on chromosome 9p (A), and 1 signal for both bio and dig probes on chromosome 3 (B).

Table 1 .Table 1 .

 11 Tumour characteristics

	Iris melanoma no#	Age years	Jakobiec group	Growth pattern	Ciliary body involved	Cell type	Angle invasion	Survival	Metastases
	1	53	9	diffuse	yes	epithelioid	yes	48	yes †
	2	63	9	diffuse	no	epithelioid	yes	112	no
	3	41	9	multiple	yes	epithelioid	yes	20	no
	4	55	9	diffuse	yes	epithelioid	yes	16	no
	5	77	6	multiple	no	spindle	yes	62	no
	6	76	6	nodular	no	spindle	no	10	no
	7	71	6	nodular	no	spindle	no	2	no
	8	42	9	multiple	no	epithelioid	yes	144	no
	9	52	8	diffuse	yes	mixed	yes	37	yes †
	10	59	8	diffuse	yes	mixed	yes	29	yes
	11	71	8	diffuse	no	mixed	yes	23	no
	12	46	7	nodular	no	spindle	no	96	no
	13	49	8	nodular	no	mixed	no	157	no
	14	62	7	multiple	no	spindle	yes	62	no
	15	69	9	diffuse	no	epithelioid	yes	5	no
	16	66	8	diffuse	yes	mixed	yes	35	yes †
	17	56	8	multiple	yes	mixed	yes	86	no
	18	69	7	nodular	no	spindle	no	24	no
	19	83	9	multiple	yes	epithelioid	yes	157	no
	20	67	7	nodular	no	spindle	no	204	no
	† Death due to metastases

Table 2 .

 2 Table 2. FISH results

	Iris melanoma					Copy number per chromosomal region
	no#	1p36 3p13	cen3	3q25	6p22	6q21	8p21	cen8	8q22	9p21	10q23
	1	1	1	1	1	2	2	2	3	3	1	2
	2	2	1	1	1	2	1	3	3	3	1	2
	13	n/a	2	2	n/a	n/a	n/a	n/a	3	n/a	n/a	n/a
	14	1	1	1	n/a	n/a	n/a	n/a	n/a	3	1	2
	15	1	1	1	1	2	2	1	3	3	2	2
	16	1	1	1	n/a	2	1	n/a	n/a	3	2	2
	17	n/a	2	2	n/a	n/a	n/a	n/a	n/a	n/a	2	2
	18	2	2	2	2	2	2	2	2	2	2	2
	19	1	1	1	n/a	n/a	n/a	n/a	n/a	3	2	2
	20	2	2	2	n/a	n/a	n/a	n/a	n/a	2	2	2
	Total aberrations	6/18 9/20	9/20	6/14	0/15	2/15	4/14	9/15	9/18	7/19	0/19

Acknowledgements

The authors are grateful to L. Razzaq for helping to complete the follow-up data in Leiden.

SWOO-Flieringa Foundation and Professor Henkes Foundation


Rowida Almomani

Jaap Van Der Heijden

Yavuz Ariyurek

Yuchig Lai

Michiel Van Galen

Martijn H Breuning

Prof. Dr Johan T Den Dunnen

Experiences with array-based sequence capture; towards clinical applications

Keywords: capture array, heterogeneous disorders, sequencing

Although sequencing of a human genome gradually becomes an option, zooming in on the region of interest remains attractive and cost saving. We performed array-based sequence capture using 385K NimbleGen arrays to zoom in on the protein-coding and immediate intron flanking sequences of 112 genes, potentially involved in mental retardation and congenital malformation. Captured material was sequenced using Illumina technology. A data analysis pipeline was built which detects sequence variants, positions them in relation to the gene, checks for presence in databases (e.g. dbSNP) and predicts the potential consequences at the level of RNA splicing and protein translation.

In the samples analyzed all known variants were reliably detected, incl. pathogenic variants from control cases and SNPs derived from array experiments. Although overall coverage varied considerably, it was reproducible per region and facilitated the detection of large deletions and duplications (CNVs), incl. a partial deletion in the B3GALTL gene from a patient sample. For ultimate diagnostic application overall results need to be improved. Future arrays should contain probes from both DNA strands and to obtain a more even coverage one could add fewer probes from densely and more probes from sparsely covered regions.

Introduction

For many years, the amplification of target sequences by PCR followed by Sanger sequencing has been the gold standard for screening of variants in terms of both read length and accuracy of sequencing. [START_REF] Sanger | DNA sequencing with chain-terminating inhibitors[END_REF] However, when it comes to conditions with highly heterogeneous etiology, a large number of different genes need to be screened for mutations. In such cases, gathering information becomes laborious, expensive and time consuming. There are many examples of diseases that can be caused by mutations in many different genes, including mental retardation (MR), [START_REF] Chelly | Genetics and pathophysiology of mental retardation[END_REF] Charcot-Marie-Tooth disease, [START_REF] Szigeti | Charcot-Marie-Tooth disease[END_REF] cardiomyopathy, [START_REF] Paul | Genes causing inherited forms of cardiomyopathies[END_REF] retinitis pigmentosa, 5 autism, [START_REF] Muhle | The genetics of autism[END_REF] hearing loss [START_REF] Hilgert | Forty-six genes causing nonsyndromic hearing impairment: which ones should be analyzed in DNA diagnostics?[END_REF] and congenital disorders of glycosylation. [START_REF] Freeze | Genetic defects in the human glycome[END_REF] Extensive re-sequencing of many disease-associated genes is required in order to explore, at the sequence and structural level, the genomic variation that might be involved in causing such diseases.

Several next generation sequencing (NGS) platforms are now available and they have allowed the sequencing and analysis of large numbers of genes in one experiment, [START_REF] Bonetta | Genome sequencing in the fast lane[END_REF][START_REF] Von Bubnoff | Next-generation sequencing: the race is on[END_REF][START_REF] Schuster | Next-generation sequencing transforms today's biology[END_REF] and are able to generate a massive amount of sequence data and have considerably reduced the cost of DNA sequencing. [START_REF] Shendure | Next-generation DNA sequencing[END_REF] However, although NGS platforms have enormously increased throughput and have permitted whole-genome sequencing, high cost still prevents routine whole human genome re-sequencing projects. Therefore, zooming in on the region of interest is an attractive option. In addition it circumvents the problem of identifying variants in genes for which the analyses was not intended (with associated ethical problems).

Microarray-based genomic selection combined with massively parallel high throughput sequencing is the method of choice to analyze large numbers of genes in a more comprehensive and cost effective way. [START_REF] Albert | Direct selection of human genomic loci by microarray hybridization[END_REF][START_REF] Okou | Microarraybased genomic selection for high-throughput resequencing[END_REF][START_REF] Hodges | Genome-wide in situ exon capture for selective resequencing[END_REF] We have used custom high-density microarrays (NimbleGen) for the enrichment of 112 distinct genes potentially involved in mental retardation and congenital malformation, followed by sequencing on the Illumina Genome Analyzer I platform.

The first aim of our study was to apply and validate the array-based enrichment method as an efficient and convenient strategy to capture any desired portion of the human genome. The second aim was to accelerate the detection of sequence and copy number variations (CNV) in the selected candidate genes with lower costs, especially for the genes that are potentially involved in MR.

Material and methods

Sample selection and validation

Six DNA samples were used in this study, including two controls containing known pathogenic variants, sample S-2 contains a known MECP2 (OMIM 300005) pathogenic point mutation (c.538C>T). The second sample, patient S-6, carry a large deletion spanning exons 8 to 15 in one allele and a splice site mutation (c.660+1G) at the other allele of the B3GALTL (OMIM 610308) gene.

The other four DNA samples were from patients with mental retardation with an unknown cause. Single nucleotide polymorphism (SNP) array data were available for two samples: S-7 with 250K Nsp Affymetrix, and S-5 with 317K Illumina data. We used these data to validate the sequences obtained after capture-array and Illumina sequencing.

Causative large deletions and duplications had been previously excluded by SNP-array testing in S-3, S-5, S-7, S-8.

Exon array design

Microarrays with 385K probe capacity (NimbleGen) were used to capture all exons, splice site and immediately adjacent intron sequence of 112 human genes. Based on searches in OMIM and literature we selected 112 human genes known to cause MR either as part of a known syndrome or in isolation Sup. Table 1. Primary sequence data from all exons was extracted from NCBI's genome (Build 36). Microarrays were designed by NimbleGen with long oligonucleotide probes (54-99 nucleotides) which span each target region, overlapped and shifted on average of seven bases. [START_REF] Albert | Direct selection of human genomic loci by microarray hybridization[END_REF] The oligonucleotides were designed to achieve isothermal hybridization across the arrays capturing one strand only.

All highly repetitive regions were excluded from the probe selection in order to avoid non-specific capturing of genomic regions. Using all criteria listed, for 2% of the target sequences no capture probe could be designed (note that theoretically these sequences can be covered partly through capture from directly flanking unique sequences). Four of the arrays were reused at least twice.

Genomic DNA library preparation and target capture

The methods used for target capture, enrichments and elution followed previously described protocols with slight modifications (Roche/NimbleGen). [START_REF] Nimblegen | NimbleGen services user's guides: sequence capture service[END_REF] Genomic DNA (20-10μg) was fragmented using a nebulizer or Bioruptor according to instructions from the manufacturer to yield fragments from 250-1000 bp (nebulization) or 250-600 bp (Bioruptor). Adapter oligonucleotides from Illumina (single reads) were ligated to the ends. After the ligation was completed, successful adapter ligation was confirmed by PCR. The DNA-adapter ligated fragments were then hybridized to the sequence capture microarray for 65 hours. After hybridization and washing, the DNA fragments bound to the array were eluted, using 300 μl of the elution buffer (Qiagen) on each array. A gasket (Agilent) was applied and placed on the thermal elution device (home-made) for 20 min at 95 ºC. We repeated this process once by adding 200 μl elution buffer (Qiagen). DNA from each eluted sample was enriched by 18-cycle PCR using a high fidelity polymerase and a single primer pair corresponding to the Illumina adapters ligated earlier.

Check enrichments by qPCR

To verify successful hybridization capture we performed qPCR (quantitative PCR) on DNA samples (S-2, S-3, S-5, S-7, S-6, S-8) before and after array enrichment. The primers amplified five loci from MBL2, DMD, and BRCA1 (100bp) as negative controls (no capture probes on the array) and four loci from MECP2, CREBBP and NSD1 genes, as positive controls (capture probes on the array) Sup. Table 2. All primers for qPCR were designed using Primer 3 (http://frodo.wi.mit.edu/). The qPCR assays were done in triplicate in the Lightcycler using 384-well plates (Roche) in 10 μl total volume: 5 μl of 2x SYBR Green master Rox (Roche), 0.25 μl of each primer (10 pmol/μl), 2 μl of DNA template and 2.5 μl of ultra-pure water. The thermo-cycling protocol was done as follows: 10 min at 95 ºC, 45 cycles of 10 s at 95 ºC, 30 s at 60 ºC, 20 s at 72 ºC, and 5 min at 72 ºC, followed by melting curve analysis in order to determine the specific and non-specific amplified products and other artifacts that might interfere with the CP values. To calculate the relative fold-enrichment of the targeted regions, we compared amplification of the positive versus negative controls. The relative fold-enrichment R was calculated using the values of Δ CP (i.e. the difference between average CP of non-captured and average CP of captured samples) according to: R = E N , Where E is the efficiency of the qPCR assay for a particular amplicon and N = Δ CP (crossing point).

DNA sequencing

The eluted-enriched DNA fragments were sequenced using the Illumina GAI platform at the Leiden Genome Technology Center (LGTC). Single-end sequencing of 36 or 50 nucleotides was performed following the instructions of the manufacturer.

Reads mapping and data analysis

Sequence read mapping was done by ELAND and ELAND-extended programs, which were a part of the Illumina GAI data analysis package. Only reads of high quality scores were mapped to the human reference genome (NCBI, BUILD 36.2), allowing up to two mismatches. We created different Perl scripts to extract and process data from the ELAND files. Coverage was calculated at the target level (gene-exons), nucleotide level and per probe region. SNP calling was performed by searching for nucleotides discordant with the reference genome with base call quality score of 30 (99.9% base call accuracy), a read depth of 8 or greater and the variant allele larger than 30% of the total coverage. Then all variants were checked for their presence in known databases, e.g. dbSNP. Perl scripts were designed to predict the potential consequences at the level of RNA splicing and protein translation based on Ensemble v.51. Furthermore, we designed a Perl script to facilitate detection of small deletions/insertions (up to three nucleotides). All Perl scripts are available upon request..

Sanger Sequencing

Twenty-one variants detected by Illumina GAI analyzer were selected and confirmed by Sanger sequencing using the standard Sanger sequencing protocol at the Leiden Genome Technology Center (LGTC). The primer sequences (with M13 tail) used are shown in the Sup. Table 3.

Results

The methodology used starts with fragmentation of the genomic DNA. Linker and primer addition can then be performed either prior to or after array-capture target enrichment. To facilitate limited amplification of the expected low yield array elution we decided to perform full Illumina sample preparation prior to array capture. Initially experiments were performed using 20μg genomic DNA, later we reduced this to 10μg. We used qPCR, comparing targeted (four positive controls) and non-targeted regions (five negative controls), to check successful array enrichment as well as to estimate the foldenrichment obtained, Sup. Tables 4 and5 for examples. Since enrichment varies significantly from locus to locus, we tested multiple loci to get an accurate estimate.

Samples where qPCR did not indicate clear enrichment (> 100x) were discarded. The ultimate enrichments achieved varied from experiment to experiment with a tendency to increase over time, indicating that lab-experience is an important aspect of the array capture technology. Since the fold-enrichments determined by qPCR correlate positively with the average sequence depth obtained, we conclude that qPCR provides an effective and cost-saving check for successful enrichment, examples are listed in Sup. Tables 4 and5.

Sequence data

The custom arrays used contained 112 different human genes that are known to be or potentially involved in mental retardation (MR) and congenital malformation. Samples were run on one channel of the Illumina GAI. For sequence analysis we used only those QC filtered reads that map back uniquely to the reference sequence (M0) or with 1 or 2 mismatches (M1, M2) (Fig 1). Using these settings 85-92% of the targeted nucleotides were covered by at least eight reads (Table 1) and 94-98% by at least one read (note that for 2% of the targeted sequences no probe could be designed, see M&M). Effectively, this means that for 78% of the targeted sequences on the array coverage was sufficient (>20x) to detect any variants that were present. Two of the samples had been previously analyzed using SNP arrays. The region selected using the capture array included 67 different SNPs that had been present on the SNP-arrays. We observed a perfect agreement (100%) between array-based SNP calls and those obtained using NGS (67/67 variants) Sup. Table 6.

To determine our ability to detect pathogenic mutations, we included one sample from a female patient (S-2) harboring a dominant pathogenic point mutation in the MECP2 gene, (c.538C>T) on the X-chromosome. Our results clearly detected the change in the heterozygous state Sup. Table 7. Similarly we detected a homozygous change in the B3GALTL gene in a Peter's Plus patient (c.660+1G>A, Sup. Table 7, see below). Next, we selected 21 variants detected in samples S-2, S-3, S-5, S-7, and S-8 and checked these by traditional Sanger sequencing. We were able to confirm 21 of the 21 variants, including their status being homozygous or heterozygous Sup. Table 7. The analysis of the variants found in all 112 genes of the patients did not reveal a clear cause of their mental retardation Sup.Table 8 and9.

Copy number variations (CNV)

Changes that cannot be easily detected using the sequence itself include deletions and duplications (CNVs). However, such variants can be expected to give quantitative changes in coverage. To see whether overall coverage can be used to detect quantitative changes, we first analysed the 39 genes located on the X-chromosome. Indeed, when coverage was normalized using autosomal genes (Fig. 2A), samples from female showed a clearly higher X-chromosome coverage compared to male samples (Fig. 2B).

Furthermore, as expected, the gene on the Y-chromosome (NLGN4Y) gave no coverage in the female sample (Fig. 2B). To determine the sensitivity of our method for detecting smaller CNVs, we carefully analysed a sample from a compound heterozygous patient (S-6) carrying a partial deletion (exons 8-15) and a splice site mutation (c.660+1G>A, intron 8) in the B3GALTL gene. The splice site mutation was evident as no wild-type sequence was present. The presence of a deletion emerged since, as compared to other samples, we observed a significantly lower average coverage for the B3GALTL gene (53x versus 155x, 150x, 140x) (Fig. 2C). In addition, while the splice site mutation in exon 8 was detected in the "homozygous" state (like all nine variants downstream), we observed variants in the first exons (1-7) also in heterozygous state Sup. Table 10. These data show that not only have we obtained an excellent specificity of the capture process but that we have also been able to distinguish between male and female samples.

Discussion

Array-based genomic selection offers several advantages for large-scale targeted DNA isolation over other approaches like PCR-based methods (long range PCR or multiplexed short PCR), [START_REF] Edwards | Multiplex PCR: advantages, development, and applications[END_REF][START_REF] Markoulatos | Multiplex polymerase chain reaction: a practical approach[END_REF][START_REF] Cutler | High-throughput variation detection and genotyping using microarrays[END_REF] selector technology [START_REF] Dahl | Multiplex amplification enabled by selective circularization of large sets of genomic DNA fragments[END_REF][START_REF] Dahl | Multigene amplification and massively parallel sequencing for cancer mutation discovery[END_REF] and BACs technology. [START_REF] Bashiardes | Direct genomic selection[END_REF] PCR-based methods become laborious, time consuming and costly if hundreds to thousands of regions (exons) need to be amplified, especially if all the sequences are required. Furthermore, when PCRs are multiplexed it becomes difficult to check successful amplification per fragment, the chance of obtaining artifacts increases and equimolar loading before sequencing becomes very difficult. New approaches for massive individual PCR have been introduced recently [START_REF] Tewhey | Microdroplet-based PCR enrichment for large-scale targeted sequencing[END_REF] but experiences with these are still limiting. Selector technology [START_REF] Dahl | Multiplex amplification enabled by selective circularization of large sets of genomic DNA fragments[END_REF][START_REF] Dahl | Multigene amplification and massively parallel sequencing for cancer mutation discovery[END_REF] seems attractive but it largely depends on proper in-house probe design and experience thus far is very limited. Successful genomic selection using BACs has been demonstrated but has several limitations. Since a BAC is the unit of selection, multiple BACs are required to isolate discontinuous regions of interest.

In this study we have tested array-based sequence capture to determine the sequence of 112 genes potentially involved in mental retardation. We show that array-based sequence capture technology is an efficient, quick and reliable method for the parallel sequencing of a range of genes of interest. Known variants (array-based calls) for 67 SNPs matched perfectly with those obtained using NGS Sup. Table 6. Two positive controls with known pathogenic changes in the MECP2 gene (sample S-2) and B3GALTL gene (sample S-6) were readily detected. In addition, 21/21 selected variants found in the 5 samples analyzed could be confirmed using Sanger sequencing Sup. Table 7. Sequence coverage of the nucleotide of interest is critical for reliably detecting sequence changes. If coverage is too low both false positives (caused by sequence errors) and false negatives (if only one allele from a heterozygous sample is observed) will occur.

The coverage we obtained differs significantly not only between targeted genomic regions (genes) but also between different samples Sup. Table 1, Fig. 2A. Since the overall methodology is rather complex, particularly the collection of the hybridized array-enriched DNA sequences, the difference between samples is most probably influenced by technical factors such as variations in the hybridization, washing conditions, and potentially re-using the capture array. Furthermore coverage is influenced by array design, including probe sequence (melting temperature, GC content), probe density, and spacing Sup. Table 1. Our data show that AT-rich regions (>55%), regions with an overall low probe density (<3) and small exons (on average 90bp) yield a low coverage, which also varies significantly between experiments. For a second-generation capturearray the results obtained could be used to change probe density, i.e. decreased in wellcovered and increased in low-covered regions.

Our data shows that longer reads (50 bp) improves accuracy and selectivity of read mapping to the reference genome which influenced the SNP calling by having less false positives and slightly better coverage.

Since CNVs (deletions / duplications) are a significant cause in the etiology of mental retardation [START_REF] Shaw-Smith | Microarray based comparative genomic hybridisation (array-CGH) detects submicroscopic chromosomal deletions and duplications in patients with learning disability/mental retardation and dysmorphic features[END_REF] we tested the feasibility of detecting large CNVs using array capture and NGS. Our results indicate that, if coverage is sufficiently high, array-capture can also be used to detect such quantitative changes. Our array contained one gene from the Ychromosome which gave no coverage in females (Fig. 2B), while the 39 X-linked genes when compared to the 69 autosomal genes yielded overall 50% lower coverage in the male samples (Fig. 2B). Another example derives from a sample containing a partial B3GALTL gene deletion on one allele (exons 8-15) and a splice site mutation on the other allele (c.660+1G>A). Although coverage over the entire gene seems reduced (experimental variation / coincidence), coverage for the second half of the gene clearly drops below that of normal (Fig. 2D). An algorithm for detecting local deviations from the average coverage is currently under development.

Regarding probe design (performed by NimbleGen) it should be noted that all array-probes are from one strand (coding DNA strand) and thus DNA molecules from only the non-coding strand are captured. This has several consequences. First, the sequence obtained is from one strand only while for diagnostic applications quality assurance requires that sequences are obtained in forward and reverse orientation.

Sequencing this one strand in both directions is partly fooling one self. Second, we observed that the sequences obtained relative to the array-probes extend in a 5' but not in a 3' direction. The most probable cause for the latter is steric hindrance during array hybridization, preventing non-hybridizing tails at the surface side of the array. When capture probes are attached with their 3' ends, this has consequences for probe design at the edges of the targeted regions; on the 5' side coverage will be significantly better than on the 3' side. Both effects could be overcome simply by reversing the probe sequence of every other nucleotide on the array. Theoretically this would also mean that the overall yield of enriched DNA would double as both strands from the sample will be captured.

To save costs, we have re-used the arrays up to three times by hybridizing different samples. The danger of this approach is of course contamination if hybridized DNA from a previous experiment is not eluted completely. Indeed, in some experiments we observed low-level contamination e.g. through heterozygous calls from Xchromosome sequences in male samples. It should be noted however that crosscontamination can be easily controlled when samples containing differently tagged linkers are used in subsequent experiments.

Using the current design, low coverage was obtained mainly at the edges of the regions targeted, especially the 3' side (see above), i.e. directly gene flanking or intronic regions. Although coverage varied widely, 78% of all regions targeted and present on the array were covered effectively by the sequence obtained. Note that there is a clear correlation between fragment size of the genomic DNA used and the coverage, the larger the fragment size used the lower the target coverage achieved, since more flanking DNA is captured. Especially for array-based capture, due to the steric hindrance described, this effect will be significant near the array-attached end of a probe targeted region. Assuming that second-generation capture-arrays will be more effective (i.e. complete and with even coverage) and sequence power will improve further, it should soon be possible to sequence-tag, mix and simultaneously analyze different samples in one experiment, giving a significant cost reduction.

Recently in-solution capture was presented as an alternative to array-based capture. [START_REF] Gnirke | Solution hybrid selection with ultra-long oligonucleotides for massively parallel targeted sequencing[END_REF] Besides advantages of simplicity, a reduced workload and a potential for automation, when attempted, in-solution capture will not show the effect of steric hindrance we observed. However, capturing both strands would be complicated by the fact that capture probes will hybridize with each other. Initial experiences in our lab with in-solution capture were successful and for future projects we will change to this approach.

Overall we conclude that array-based sequence capture followed by next generation sequencing offers a versatile tool for successfully selecting sequences of interest from a total human genome. The approach will be especially helpful in speeding up the identification of the pathogenic mutation(s) in diseases where the genomic region to be scanned is large. Our results indicate that the methodology can still be improved, in particular with respect to probe design obtaining a more even coverage of the targeted regions. Based on initial experiences and publications, we expect that array-capture will quickly be replaced by in-solution capture. Ultimately, the cost of this approach is determined by the minimal coverage, which in turn determines the sensitivity required for the detection of potential sequence variants. shows average coverage of 69 autosomal genes from four different samples. (B) shows average coverage of 39 genes located on X and one gene (NLGN4Y) located on the Ychromosome, a female sample exhibited an absence of hybridization on the captured array, with 0 coverage in the regions corresponding to the NLGN4Y. The female sample shows a higher average coverage per gene for all genes located on X-chromosome compared to male samples. (C) Lower average coverage of B3GALTL gene in a male patient sample with a known large deletion compared with three wild type male samples. (D) Average coverage per nucleotide for the second half (exons 8-15) of B3GALTL gene, the patient sample shows lower average coverage for this region compared with wild type. del = deletion, wt = wild type.

Table and

Table 1 :

 1 Sequence summary results of the different array-capture experiments performed. QC = quality control, MM# reads = number of reads with # mismatches to the reference sequence, F = female, M =male.

Figure 1 .

 1 Figure 1. Detection of sequence variants. 32 nucleotide NGS reads (top, sequence

Figure 2 .

 2 Figure 2. Average coverage obtained for different genes in four different samples. (A)

ACKNOWLEDGEMENTS

We would like to thank the Leiden Genome Technology Center (LGTC), in particular Sophie Greve-Onderwater, Matthew Hestand and Rolf Vossen for their expert technical assistance. Antoinette Gijsbers for sharing the SNP data, and Kamlesh Madan for critical reading of the manuscript. The research leading to these results has received funding from the European Community's Seventh Framework Program (FP7/2007-2013) under grant agreements 223026 (NMD-chip) and 223143 (the TechGene).