Principal component analysis (PCA) was applied for evaluating the selectivity of 22 GC columns for which complete retention data were available for the 136 tetra- to octa-chlorinated dibenzo-p-dioxins (PCDDs) and dibenzofurans (PCDFs). Because the hepta- and octa-homologues are easy to separate the PCA was focused on the 128 tetra- to hexa-CDD/Fs. The analysis showed that 21 of the 22 GC columns could be subdivided into four groups with different selectivity. Group I consists of columns with non-polar thermally stable phases (Restek 5Sil MS and Dioxin 2, SGE BPX-DXN, Supelco Equity-5, and Agilent DB-1, DB-5, DB-5ms, VF-5ms, VF-Xms and DB-XLB). Group II includes ionic liquid columns (Supelco SLB-IL61, SLB-IL111 and SLB-IL76) with very high polarity. Group III includes columns with high-percentage phenyl and cyanopropyl phases (Agilent DB-17 and DB-225, Quadrex CPS-1, Supelco SP-2331, and Agilent CP-Sil 88), and Group IV columns with shape selectivity (Dionex SB-Smectic and Restek LC-50, Supelco βDEXcst, Agilent VF-Xms and DB-XLB). Thus, two columns appeared in both Group I and IV (Agilent VF-Xms and DB-XLB). The selectivity of the other column, Agilent DB-210, differs from those of these four groups. Partial least squares (PLS) regression was used to correlate the retention times of the tetra- to hexa-CDD/Fs on the 22 stationary phases with a set of physicochemical and structural descriptors to identify parameters that significantly influence the solute-stationary phase interactions. The most influential physicochemical parameters for the interaction were associated with molecular size (as reflects in the total energy, electron energy, core-core repulsion and standard entropy), solubility (aqueous solubility and n-octanol/water partition coefficient), charge distribution (molecular polarizability and dipolar moment), and reactivity (relative Gibbs free energy); and the most influential structural descriptors were related to these parameters, in particular, size and dipolar moment. Finally, the PCA and PLS analyses were complemented with linear regression analysis to identify the most orthogonal column combinations, which could be used in comprehensive two-dimensional gas chromatography (GC×GC) to enhance PCDD/F separation and congener profiling.